Variational Amplitude Amplification for Solving QUBO Problems

1. Introduction

Amplitude amplification is a quantum algorithm strategy that is capable of circumventing one of quantum computing’s most difficult challenges: probabilistic measurements. Originally proposed by Grover in 1996 [1], and later shown to be optimal [2,3], the combination of his oracle $U_{\mathrm{G}}$ and `diffusion’ $U_{\mathrm{s}}$ operators is able to drive a quantum system to a superposition state where one (or multiple) basis state(s) has nearly 100% probability of being measured. Since then, many researchers have contributed to the study of $U_{\mathrm{G}}$ and $U_{\mathrm{s}}$ [4,5,6,7,8,9], seeking to better understand how the fundamental nature of amplitude amplification is dependent on these two operators. Similarly, the aim of this study is to further extend the capabilities of amplitude amplification as a means for solving combinatorial optimization problems using gate-based quantum computers.

The results of this paper are a continuation of our previous work [10], in which we demonstrated an oracle design which was capable of encoding and solving a weighted directed graph problem. The motivation for this oracle was to address a common criticism of $U_{\mathrm{G}}$[11,12,13,14,15], namely that the circuit construction of oracles too often hardcodes the solution it aims to find, negating the use of quantum entirely. Similar to other recent studies [16,17,18,19], we showed that this problem can be solved at the circuit depth level by avoiding gates such as control-Z for constructing the oracle, and instead using phase and control-phase gates (P($\theta$) and CP($\theta$)). However, simply changing the phase produced from $U_{\mathrm{G}}$ to something other than $\pi$ is not enough [20,21,22,23,24,25]. Our oracle construction applies phases to not only a desired marked state(s), but $all$ states in the full $2^N$ Hilbert Space. The phase each basis state receives is proportional to the solutions of a weighted combinatorial optimization problem, for which the diffusion operator $U_{\mathrm{s}}$ can be used to boost the probability of measuring states that correspond to optimal solutions.

The consequence of using an oracle operation that applies phases to every basis state is an interesting double-edged sword. As we show in sections II. - IV., and later in section VII., the use of phase gates allows for amplitude amplification to encode a broad scope of combinatorial optimization problems into oracles, which we call `cost oracles’ $U_{\mathrm{c}}$. In particular, we demonstrate the robustness of amplitude amplification for solving these kinds of optimization problems with asymmetry and randomness [26,27,28]. However, the tradeoff for solving more complex problems is twofold. Firstly, in contrast to Grover’s oracle, using $U_{\mathrm{c}}$ is only able to achieve peak measurement probabilities up to 70-90%. In section VI. we show that these probabilities are still high enough for quantum to reliably find optimal solutions, which notably are achieved using the same O( $\frac{\pi }{4}\sqrt{N/M}$ ) iterations as standard Grover’s [1,2,3].

The second, more challenging tradeoff when using $U_{\mathrm{c}}$ is that the success of amplitude amplification is largely dependant on the correct choice of a single free parameter $p_{\mathrm{s}}$[10]. This scalar parameter is multiplied into every phase gate for the construction of $U_{\mathrm{c}}$ (P($\theta ·p_{\mathrm{s}}$) and CP($\theta ·p_{\mathrm{s}}$)), and is responsible for transforming the numeric scale of a given optimization problem to values which form a range of approximately $2\pi$. This in turn is what allows for reflections about the average amplitude via $U_{\mathrm{s}}$ to iteratively drive the probability of desired solution states up to 70-90%. The significance of $p_{\mathrm{s}}$, and the challenges in determining it experimentally, are a major motivation for this study. In particular, the results of section V. demonstrate that there is a range of $p_{\mathrm{s}}$ values for which many optimal solutions can be made to become highly probable. Additionally, our simulations show that there is an observed correlation between the numerical cost function value of these solutions and the $p_{\mathrm{s}}$ values where they achieve peak probabilities. This underlying correlation supports the idea of using amplitude amplification for a variational model of hybrid quantum-classical computing, which is the core finding of this study.

1.1. Layout

The layout of this study is as follows. Section II. begins with the mathematical formalism for the optimization problem we will seek to solve using amplitude amplification. Sections III. & IV. discuss the construction of the problem as a quantum circuit, the varying degrees of success one can expect from optimization problems generated using random numbers, and the conditions for which these successes can be experimentally realized. In section V. we explore the role of $p_{\mathrm{s}}$ from a heuristic perspective, whereby we demonstrate that many near optimal solutions are capable of reaching significant probabilities of measurement. Section VI. is a primarily speculative discussion, theorizing how the collective results of section V. can be coalesced into a hybrid quantum-classical variational algorithm. And finally, section VII. completes the study with additional optimization problems that can be constructed as oracles and solved using amplitude amplification.

2. QUBO Definitions

We begin by outlining the optimization problem which will serve as the focus for this study: QUBO (quadratic unconstrained binary optimization). The QUBO problem has many connections to important fields of computer science [29,30,31,32,33], making it relevant for demonstrating quantum’s potential for obtaining solutions. To date, the two most successful quantum approaches to solving QUBOs are annealing [34,35,36,37] and QAOA [38,39,40,41], with a lot of interest in comparing the two [42,43,44]. Shown below in Equation (1) is the QUBO cost function C(X) which we shall seek to solve using our quantum algorithm.

$$\begin{array}{c}\hfill \mathrm{C}\left(\mathrm{X}\right)=\sum _i^N{W}_i{x}_i+\sum _{\{i,j\}\in \mathbb{S}}{w}_{ij}{x}_i{x}_j\end{array}$$

(1)

The function C(X) evaluates a given binary string X of length N, composed of individual binary variables $x_i$. Together, the total number of unique solutions to each QUBO is $2^N$, which is also the number of quantum states producible from N qubits. Throughout this study we will use subscripts X${}_i$ and C(X${}_i$) when referring to individual solutions, and C(X) when discussing a cost function more generally.

As shown in Equation (1), a QUBO is defined by two separate summations of weighted values. The first summation evaluates weights $W_i$ associated with each individual binary variable, while the second summation accounts for pairs of variables which share a weighted connection $w_{ij}$. In this study we adopt the typical interpretation of QUBOs as graph problems, whereby each binary variable $x_i$ represents a node. We can then define the connectivity of a QUBO graph using the set $\mathbb{S}$, which itself is a collection of sets that describe each pair of nodes $x_i$ and $x_j$ that share a connection. See Figure 1 below for an example.

The interest of this study is to use a quantum algorithm to find either X${}_{\min }$ or X${}_{\max }$, which are the solutions which minimize / maximize the cost function C(X) respectively. For all QUBOs analyzed in the coming sections, the weight values $W_i$ and $w_{ij}$ are restricted to integers, randomly selected from a uniform distribution as shown below in Equations (2) and (3).

$$\begin{array}{ccc}\hfill W_i,w_{ij}& \in & \mathbb{Z}\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill W_i,w_{ij}& \in & [-100,100]\hfill \end{array}$$

(3)

In section V. we discuss the consequences of choosing weight values in this manner and its advantage for quantum. However, nearly all of the results shown throughout this study are applicable to the continuous cases for $W_i$ and $w_{ij}$ as well, with the one exception being the results of section V.D.

2.1. Linear QUBO

The cost function given in Equation (1) is applicable to any graph structure $\mathbb{S}$, so long as every node and edge is assigned a weight. For this study we will focus on one specific $\mathbb{S}$, which we refer to as a `linear QUBO.’ The connectivity of these graphs is as follows:

$$\begin{array}{ccc}\hfill \mathbb{S}& =& \left\{\right\{n,n+1\left\}\right|1\le n\le N-1\}\hfill \end{array}$$

(4)

As the name suggests, linear QUBOs are graphs for which every node has connectivity with exactly two neighboring nodes, except for the first and final nodes. The motivation for studying QUBOs of this nature is their efficient realizability as quantum circuits, given in the next section.

3. Amplitude Amplification

The quantum strategy for finding optimal solutions to C(X) investigated in this study is amplitude amplification [4,5,6,7,8,9], which is the generalization of Grover’s algorithm [1]. The full algorithm is shown below in Algorithm 1, which notably is almost identical to Grover’s algorithm except for the replacement of Grover’s oracle $U_{\mathrm{G}}$ with our cost oracle $U_{\mathrm{c}}$ .

Algorithm 1 Amplitude Amplification Algorithm

Initialize Qubits: $|\mathsf{\Psi }\rangle ={|0\rangle }^{\otimes N}$  Prepare Equal Superposition: $H^{\otimes N}|\mathsf{\Psi }\rangle =|s\rangle$  for$k\approx \frac{\pi }{4}\sqrt{2^N}$do       Apply $U_{\mathrm{c}}|\mathsf{\Psi }\rangle$ (Cost Oracle)      Apply $U_{\mathrm{s}}|\mathsf{\Psi }\rangle$ (Diffusion)  end for  Measure

By interchanging different oracle operations into the Algorithm 1, various problem types can be solved using ampltitude amplification. For example, Grover’s original oracle solves an unstructured search, whereas here we are interested in optimal solutions to a cost function. Later in section VII. we discuss further oracle adaptations and the problems they solve. For all oracles, we use the standard diffusion operator $U_{\mathrm{s}}$, given below in Equation (5).

$$\begin{array}{c}\hfill U_{\mathrm{s}}=2|s\rangle \langle s|-\mathbb{I}\end{array}$$

(5)

This operation achieves a reflection about the average amplitude, whereby every basis state in $|\mathsf{\Psi }\rangle$ is reflected around their collective mean in the complex plane. This operation causes states’ distance from the origin to increase or decrease based on their location relative to the mean, which in turn determines their probability of measurement. Therefore, a successful amplitude amplification is able to drive the desired basis state(s) as far from the origin as possible, up to a maximum distance of 1 (measurement probability of 100%).

3.1. Solution Space Distribution

A prerequisite for the success of amplitude amplification as demonstrated in this study is an optimization problem’s underlying solution space distribution; that is, the manner in which all possible solutions to the problem are distributed with respect to one another. For QUBOs, these are the $2^N$ possible C(X${}_i$) cost function values. Shown below in Figure 2 is a histogram of one such solution space distribution, for the case of a length 20 linear QUBO according to Equations (1) - (4). The x-axis represents all possible cost function evaluations, and the y-axis is the corresponding number of unique X${}_i$ solutions that result in the same C(X${}_i$) value.

Depicted in Figure 2 are all $2^{20}$ possible solutions to an example linear QUBO. Because this QUBO was generated from randomized weights, the combination of the Law of Large Numbers [45] and Central Limit Theorem [46] predicts that its underlying solution space should be approximately gaussian [47] in shape, given by Equation (6).

$$\begin{array}{c}\hfill \mathrm{G}\left(x\right)=\alpha {\mathrm{e}}^{\frac{{(x-\mu )}^2}{2{\sigma }^2}}\end{array}$$

(6)

Indeed, the histogram shown is approximately gaussian, but importantly it has imperfections resulting from the randomized weights. At large enough problem sizes (around $N\ge 20$), these imperfections have minimal impact on a problem’s aptitude for amplitude amplification, which was a result from our previous study [10]. Similarly, another recent study [18] demonstrated that in addition to symmetric gaussians, solution space distributions for both skewed gaussians and exponential profiles also lead to successful amplitude amplifications. The commonality between these three distribution shapes is that they all possess large clusters of solutions that are sufficiently distanced from the optimal solutions we seek to boost. This can be seen in Figure 2 as the location of X${}_{\min }$ and X${}_{\max }$ as compared to the central peak of the gaussian. When appropriately encoded as an oracle $U_{\mathrm{c}}$, these clusters serve to create a mean point in the complex plane which the optimal solution(s) use to reflect about and increase in probability.

3.2. Cost Oracle $U_{\mathrm{c}}$

In order to use Algorithm 1 for finding the optimal solution to a given cost function, we must construct a cost oracle $U_{\mathrm{c}}$ which encodes the weighted information and connectivity of the problem. In our previous study we referred to this operation as a `phase oracle’ $U_{\mathrm{P}}$[10], and similarly it has also been called a `subdivided phase oracle’ SPO [17,18] or `non-boolean oracle’ [19]. How one constructs $U_{\mathrm{c}}$ is problem specific, but the general strategy is to primarily use two quantum gates, shown below in Equations (7) and (8).

$$\begin{array}{ccc}\hfill \mathrm{P}\left(\theta \right)& =& \left[\begin{array}{cc}1& 0\\ 0& {\mathrm{e}}^{i\theta }\end{array}\right]\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill \mathrm{CP}\left(\theta \right)& =& \left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& {\mathrm{e}}^{i\theta }\end{array}\right]\hfill \end{array}$$

(8)

The single and two-qubit gates P($\theta$) and CP($\theta$) are referred to as phase gates, also known as R${}_{\mathrm{z}}$($\theta$) and CR${}_{\mathrm{z}}$($\theta$) for their effect of rotating a qubit’s state around the z-axis of the Bloch sphere. Mathematically they are capable of applying complex phases as shown below.

$$\begin{array}{ccc}\hfill \mathrm{P}\left(\theta \right)|1\rangle & =& {\mathrm{e}}^{i\theta }|1\rangle \hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill \mathrm{CP}\left(\theta \right)|11\rangle & =& {\mathrm{e}}^{i\theta }|11\rangle \hfill \end{array}$$

(10)

Applying P($\theta$) to a qubit only affects the $|1\rangle$ state, leaving $|0\rangle$ unchanged, and similarly only $|11\rangle$ for CP($\theta$). However, this is exactly what we need in order to construct C(X) from Equation (1). When evaluating a particular binary string X${}_i$ classically, only instances where the binary values $x_i$ are equal to 1 yield non-zero terms in the summations. For quantum, each binary string X${}_i$ is represented by one of the $2^N$ basis states $|{\mathrm{X}}_i\rangle$. Thus, our quantum cost oracle $U_{\mathrm{c}}$ can replicate C(X) by using P($\theta$) and CP($\theta$) to only effect basis states with qubits in the $|1\rangle$ and $|11\rangle$ states.

Shown above in Figure 3 is an example of a 4-qubit QUBO cost oracle, where the weighted values $W_i$ and $w_{ij}$ are used as the $\theta$ parameters for the various phase gates. Although incomplete, we will use this oracle circuit to demonstrate quantum’s ability to encode a cost function C(X). For example, consider the binary solution X${}_i=1101$ and the corresponding quantum basis state $|1101\rangle$. The classical evaluation of this solution is as follows:

$$\begin{array}{ccc}\hfill \mathrm{C}\left(1101\right)& =& -8+18-22-12\hfill \\ & =& -24\hfill \end{array}$$

(11)

Now let us compare this to the phase of $|1101\rangle$ after applying $U_{\mathrm{c}}$:

$$\begin{array}{ccc}\hfill U_{\mathrm{c}}|1101\rangle & =& {\mathrm{e}}^{i(-8+18-22-12)}|1101\rangle \hfill \\ & =& {\mathrm{e}}^{-24i}|1101\rangle \hfill \end{array}$$

(12)

The phase acquired in Equation (12) is equivalent to the classical evaluation shown in (11), which means that $U_{\mathrm{c}}$ is an accurate encoding of C(X). If we were to now apply $U_{\mathrm{c}}$ to the equal superposition state $|\mathrm{s}\rangle$ (step 2 in Algorithm 1), all $2^N$ basis states would receive phases equal to their cost function value. This is the advantage that quantum has to offer: simultaneously evaluating all possible solutions of a cost function through superposition.

3.3. Scaling Parameter $p_{\mathrm{s}}$

While the cost oracle shown in Figure 3 is capable of reproducing C(X), its use in algorithm 1 will not yield the optimal solution X${}_{\min }$ or X${}_{\max }$. This is because quantum phases are $2\pi$ modulo, which is problematic if the numerical scale of C(X) exceeds a range of $2\pi$. Consequently if two quantum states receive phases that differ by a multiple of $2\pi$, then they will both undergo the amplitude amplification process identically. If this happens unintentionally via $U_{\mathrm{c}}$, then our cost oracle cannot be used to minimize or maximize C(X).

In order to construct $U_{\mathrm{c}}$ such that it is usable for amplitude amplification, a scalar parameter $p_{\mathrm{s}}$ must be included in all of the phase gates. The value of $p_{\mathrm{s}}$ is problem specific, but its role is always the same: scaling the cumulative phases applied by $U_{\mathrm{c}}$ down (or up) to a range where [C(X${}_{\min }$) , C(X${}_{\max }$)] is approximately [x , x+$2\pi$]. This range does not have to be [0 , $2\pi$] exactly, so long as the phases acquired by $|{\mathrm{X}}_{\min }\rangle$ and $|{\mathrm{X}}_{\max }\rangle$ are roughly $2\pi$ different. See Figure 4 below for an example of $p_{\mathrm{s}}$ in $U_{\mathrm{c}}$’s construction.

Using the scaled oracle shown in Figure 4 above, let us now show how this new $U_{\mathrm{c}}$ acts on the basis state $|1101\rangle$ from before:

$$\begin{array}{ccc}\hfill U_{\mathrm{c}}|1101\rangle & =& {\mathrm{e}}^{i(-8·p_{\mathrm{s}}+18·p_{\mathrm{s}}-22·p_{\mathrm{s}}-12·p_{\mathrm{s}})}|1101\rangle \hfill \\ & =& {\mathrm{e}}^{i(-8+18-22-12)·p_{\mathrm{s}}}|1101\rangle \hfill \\ & =& {\mathrm{e}}^{-24i·p_{\mathrm{s}}}|1101\rangle \hfill \end{array}$$

(13)

As shown in Equation (13) above, multiplying $p_{\mathrm{s}}$ into every phase gate has the net effect of scaling the cumulative phase applied by $U_{\mathrm{c}}$: e${}^{-24i}\to$ e${}^{-24i·p_{\mathrm{s}}}$. Note that this is $not$ a global phase, which would have an additive effect on all states rather than a multiplicative one like shown above.

Finding the optimal $p_{\mathrm{s}}$ value for boosting X${}_{\min }$ or X${}_{\max }$ is non-trivial, and was a major focus of our previous study [10], as well as this one. In general, the scale of $p_{\mathrm{s}}$ needed for finding the optimal solution can be obtained using Equation (14) below, which scales the numerical range of a problem [C(X${}_{\min }$) , C(X${}_{\max }$)] to exactly [x , x+$2\pi$].

$$\begin{array}{c}\hfill p_{\mathrm{s}}=\frac{2\pi }{\mathrm{C}\left({\mathrm{X}}_{\max }\right)-\mathrm{C}\left({\mathrm{X}}_{\min }\right)}\end{array}$$

(14)

Although Equation (14) above is guaranteed to solve the $2\pi$ modulo phase problem mentioned previously, it is almost never the $p_{\mathrm{s}}$ value which can be used to find X${}_{\min }$ or X${}_{\max }$. Only in the case of a perfectly symmetric solution space distribution is Equation (14) the optimal $p_{\mathrm{s}}$ value, in which case the states $|{\mathrm{X}}_{\min }\rangle$ and $|{\mathrm{X}}_{\max }\rangle$ undergo the amplitude amplification process together. However, realistic optimization problems can be assumed to have a certain degree of randomness or asymmetry to their solution space, producing distributions more akin to Figure 7. For this reason, Equation (14) is better thought of as the starting point for finding the true optimal $p_{\mathrm{s}}$, which we discuss later in section IV.B. For now, Equation (14) is sufficient for demonstrating $p_{\mathrm{s}}$’s role in creating an average amplitude suitable for boosting $|{\mathrm{X}}_{\min }\rangle$ or $|{\mathrm{X}}_{\max }\rangle$, shown in Figure 5.

The bottom plot in Figure 5 shows $|\mathsf{\Psi }\rangle$ after the first application of $U_{\mathrm{c}}$ in Algorithm 1. Note the location of the average amplitude (red `x’), which is only made possible by the majority of quantum states which recieve phases near the center of the gaussian in the top plot. Optimal amplitude amplification occurs when the desired state for boosting is exactly $\pi$ phase different with the mean [2,3], which is very close to the situation seen in Figure 5. However, since this $U_{\mathrm{c}}$ is derived from a QUBO with randomized weights, the $p_{\mathrm{s}}$ value provided from Equation (14) does not exactly produce a $\pi$ phase difference between the optimal states (black star) and the mean amplitude (red `x’). Consequently, the state(s) which does become highly probable from amplitude amplification for this particular $p_{\mathrm{s}}$ is not $|{\mathrm{X}}_{\min }\rangle$ and $|{\mathrm{X}}_{\max }\rangle$, which will be the subject of the coming two sections.

4. Gaussian Amplitude Amplification

The amplitude space plot depicted at the bottom of Figure 5 is useful for visualizing how a gaussian solution space distribution can be used for boosting, but the full amplitude amplification process is far more complicated. This is especially true for the QUBOs of this study, which are generated with randomized weights. Consequently, all of the results which follow throughout the remainder of this study are produced from classical simulations of amplitude amplification using cost oracles derived from linear QUBOs according to Equations (1) - (4). For a deeper mathematical insight into these processes, please see [16,17,18].

4.1. Achievable Probabilities

Amplitude amplifiation is an appealing quantum algorithm because it solves one of the most fundamental problems of quantum computing: measurement probability. For example, a single marked state using Grover’s oracle with 30 qubits is capable of achieving a final probability that is only less than $100\%$ by one billionth of a percent [1]. Thus, a natural question to ask when using $U_{\mathrm{c}}$ is what kinds of probabilities can it produce for $|{\mathrm{X}}_{\min }\rangle$ or $|{\mathrm{X}}_{\max }\rangle$? To answer this we conducted a statistically study of linear QUBOs ranging from length $N=17$ to 27. For each N we generated numerous QUBOs according to Equations (1) - (4), totals given in appendix $\mathbf{A}$. We then let a classical simulator find the $p_{\mathrm{s}}$ value which maximized the probability of measuring $|{\mathrm{X}}_{\min }\rangle$ for each QUBO (and for certain cases the optimal $p_{\mathrm{s}}$ for $|{\mathrm{X}}_{\max }\rangle$ aswell). Results for each problem size are shown below in Figure 6.

$$\begin{array}{ccc}\hfill \mu & =& \frac{1}{2^N}\sum _i^{2^N}\mathrm{C}\left({\mathrm{X}}_i\right)\hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill \sigma & =& \sqrt{\frac{{\sum }_i^{2^N}{(\mathrm{C}\left({\mathrm{X}}_i\right)-\mu )}^2}{2^N}}\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill {\sigma }^{\prime }& =& \sigma ·p_{\mathrm{s}}\hfill \end{array}$$

(17)

Figure 6 tracks three noteworthy trends found across the various QUBO sizes: the average peak probability achievable for $|{\mathrm{X}}_{\min }\rangle$ (black triangle), the highest recorded probability for $|{\mathrm{X}}_{\min }\rangle$ (red star), and the average scaled standard deviation ${\sigma }^{\prime }$ (blue circle). For clarity, the derivation of ${\sigma }^{\prime }$ is given by Equations (15) - (17). This quantity is the standard deviation of a QUBO’s solution space distribution after being scaled by $p_{\mathrm{s}}$, making it a comparable metric for all QUBO sizes. In our previous study we demonstrated a result in agreement with Figure 6, which is the correlation between higher achievable probabilities for $|{\mathrm{X}}_{\min }\rangle$ (red star) and smaller scaled standard deviations ${\sigma }^{\prime }$ (blue circle) [10]. The latter is what is responsible for increasing the distance between $|{\mathrm{X}}_{\min }\rangle$ and the average amplitude like shown in Figure 5.

4.2. Solution Space Skewness

The relation between N, ${\sigma }^{\prime }$, and highest prob.($|{\mathrm{X}}_{\min }\rangle$) from Figure 6 can be summarized as follows: larger problem sizes tend to produce smaller standard deviations, which in turn lead to better probabilities produced from amplitude amplification. However, there is a very apparent disconnect between the probabilities capable of each problem size (red stars) versus the average (black triangle). To explain this, we must first introduce the quantity X${}_{\Delta }$ given in Equation (18) below.

$$\begin{array}{ccc}\hfill {\mathrm{X}}_{\Delta }& =& 2\mu -(\mathrm{C}\left({\mathrm{X}}_{\max }\right)+\mathrm{C}\left({\mathrm{X}}_{\min }\right))\hfill \end{array}$$

(18)

The quantity X${}_{\Delta }$ from Equation (18) is the difference between $\mathrm{C}\left({\mathrm{X}}_{\min }\right)$ and $\mu$ (the mean) minus the difference between $\mu$ and $\mathrm{C}\left({\mathrm{X}}_{\max }\right)$. A positive value for X${}_{\Delta }$ indicates that the mean is closer to $\mathrm{C}\left({\mathrm{X}}_{\max }\right)$, and vice versa for a negative valued X${}_{\Delta }$. In essence, it is a measure of skewness that describes the assymetry of a solution space distribution. Figure 7 shows example QUBO distributions for three cases of X${}_{\Delta }$, for $N=25$, demonstrating the impact X${}_{\Delta }$ has on the ability to boost $|{\mathrm{X}}_{\min }\rangle$ versus $|{\mathrm{X}}_{\max }\rangle$. While ${\sigma }^{\prime }$ is a strong indicator of a problem’s overall aptitude for amplitude amplification, X${}_{\Delta }$ determines whether the optimal minimum or maximum solution is boostable, and which is not. Further evidence of this can be seen in Figure 8, which shows 1000 randomly generated linear QUBOs of length $N=23$, and the peak probabilities achievable for $|{\mathrm{X}}_{\min }\rangle$ and $|{\mathrm{X}}_{\max }\rangle$ as a function of X${}_{\Delta }$.

Figure 7. Three randomly generated QUBO distributions for $N=25$, illustrating X${}_{\Delta }$ cases for largely positive (top), largely negative (middle), and near zero (bottom). In all three plots the exact X${}_{\Delta }$ value is reported, as well as the highest achievable probability for $|{\mathrm{X}}_{\min }\rangle$ and $|{\mathrm{X}}_{\max }\rangle$ (each from a different $p_{\mathrm{s}}$ value). Also shown in each plot are the values for C(X${}_{\min }$) and C(X${}_{\max }$), and their numerical distance to the mean $\mu$ (red-dashed line).

Figure 7. Three randomly generated QUBO distributions for $N=25$, illustrating X${}_{\Delta }$ cases for largely positive (top), largely negative (middle), and near zero (bottom). In all three plots the exact X${}_{\Delta }$ value is reported, as well as the highest achievable probability for $|{\mathrm{X}}_{\min }\rangle$ and $|{\mathrm{X}}_{\max }\rangle$ (each from a different $p_{\mathrm{s}}$ value). Also shown in each plot are the values for C(X${}_{\min }$) and C(X${}_{\max }$), and their numerical distance to the mean $\mu$ (red-dashed line).

If we compare the average peak probabilites for $|{\mathrm{X}}_{\min }\rangle$ from Figure 6 with the full data of QUBOs shown in Figure 8, we can see why the average peak probability is significantly lower than the highest recorded. Across the 1000 QUBOs studied, it is clear that X${}_{\Delta }=0$ is a dividing point for whether $|{\mathrm{X}}_{\min }\rangle$ or $|{\mathrm{X}}_{\max }\rangle$ is capable of reaching a significant probability of measurement through amplitude amplification. For $N=23$, the average prob.($|{\mathrm{X}}_{\min }\rangle$) reported in Figure 6 is approximately 64%. However, if instead we only consider QUBOs with X${}_{\Delta }>0$ from Figure 8, then the average peak probability for $|{\mathrm{X}}_{\min }\rangle$ is around 86%, and likewise for $|{\mathrm{X}}_{\max }\rangle$ when X${}_{\Delta }<0$.

Together, Figure 7 and Figure 8 demonstrate the significance of knowing X${}_{\Delta }$ from an experimenter’s perspective. Depending on the optimization problem of interest, it is reasonable to assume that an experimenter may be interested in finding only X${}_{\min }$ or X${}_{\max }$. But without any $\mathit{a}$$\mathit{priori}$ knowledge of a problem’s underlying solution space, specifically X${}_{\Delta }$, the experimenter may unknowingly be searching for a solution which is probabilistically near impossible to find through amplitude amplification. For example, consider the QUBO distribution illustrated in the top plot of Figure 7, and the peak probability for boosting $|{\mathrm{X}}_{\max }\rangle$: $0.16$%. Although it is ideal to have insight into a particular problem’s X${}_{\Delta }$ before using amplitude amplification, as we demonstrate in section V., information about X${}_{\Delta }$ can be inferred through measurement results.

4.3. Sampling for $p_{\mathrm{s}}$

If a particular optimization problem is suitable for amplitude amplification, then the speed of the quantum algorithm outlined in this study is determined by how quickly the optimal $p_{\mathrm{s}}$ value can be found. Here we shall show that sampling a cost function C(X) can provide reliable information for approximating $p_{\mathrm{s}}$ from Equation (14), which can then be used to begin the variational approach outlined in sections V. and VI. Importantly, the number of cost function evaluations needed is significantly less than either a classical or quantum solving speed. The strategy outlined in Equations (19) - (29) below can be used for approximating $p_{\mathrm{s}}$ when the experimenter is expecting an underlying solution space describable by a gaussian function (Equation (6)). If another type of distribution is expected, then the function used in Equation (22) could in principle be modified accordingly (for example, sinusoidal, polynomial, exponential [18]).

Suppose we sample a particular cost function C(X) M times, where $M<<2^N$. We will define the set $\mathbb{M}$ as the collection of values C(X${}_i$) obtained from these samples.

$$\begin{array}{c}\hfill \mathbb{M}=\{\mathrm{C}\left({\mathrm{X}}_1\right),\mathrm{C}\left({\mathrm{X}}_2\right),...,\mathrm{C}\left({\mathrm{X}}_M\right)\}\end{array}$$

(19)

Using these M values, we can compute an approximate mean and standard deviation.

$$\begin{array}{cc}\hfill \tilde{\mu }& =\frac{1}{M}\sum _{c\in \mathbb{M}}c\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill \tilde{\sigma }& =\sqrt{\frac{{\sum }_{c\in \mathbb{M}}{\left(c-\tilde{\mu }\right)}^2}{M}}\hfill \end{array}$$

(21)

In order to use Equation (14) for obtaining $p_{\mathrm{s}}$, we need approximations for C(X${}_{\min }$) and C(X${}_{\max }$). If we assume an underlying gaussian structure to the problem’s solution space, then we can write down the following equation to describe it:

$$\begin{array}{ccc}\hfill 2^N& =& {\int }_{-\infty }^{\infty }\tilde{\alpha }{\mathrm{e}}^{\frac{{(x-\tilde{\mu })}^2}{2{\tilde{\sigma }}^2}}dx\hfill \end{array}$$

(22)

$$\begin{array}{ccc}& =& -\tilde{\alpha }\tilde{\sigma }\sqrt{\frac{\pi }{2}}\mathrm{erf}{\left(\frac{\tilde{\mu }-x}{\sqrt{2}\tilde{\sigma }}\right)}_{-\infty }^{\infty }\hfill \end{array}$$

(23)

$$\begin{array}{ccc}& =& -\tilde{\alpha }\tilde{\sigma }\sqrt{\frac{\pi }{2}}·[-1-1]\hfill \end{array}$$

(24)

where erf() is the gaussian error function. Using Equation (24), we can rearrange terms and solve for an approximation to the height of the gaussian.

$$\begin{array}{c}\hfill \tilde{\alpha }=\frac{2^{N-1}}{\tilde{\sigma }\sqrt{\frac{\pi }{2}}}\end{array}$$

(25)

With the values $\tilde{\mu }$, $\tilde{\sigma }$, and $\tilde{\alpha }$ obtained from sampling, we can now approximate C(X${}_{\min }$) and C(X${}_{\max }$) using Equation (26) below.

$$\begin{array}{c}\hfill \tilde{\mathrm{G}}\left(x\right)=\tilde{\alpha }{\mathrm{e}}^{\frac{{(x-\tilde{\mu })}^2}{2{\tilde{\sigma }}^2}}=1\end{array}$$

(26)

Solving for x yields the following two values:

$$\begin{array}{c}\hfill x_{\pm }=\tilde{\mu }\pm \tilde{\sigma }\sqrt{-2\ln \left(\frac{1}{\tilde{\alpha }}\right)}\end{array}$$

(27)

which can be expressed in terms of the two quantities originally derived from sampling:

$$\begin{array}{c}\hfill x_{\pm }=\tilde{\mu }\pm \tilde{\sigma }\sqrt{-2\ln \left(\frac{\tilde{\sigma }\sqrt{\pi /2}}{2^{N-1}}\right)}\end{array}$$

(28)

And finally, the solutions $x_{\pm }$ can be used to obtain $p_{\mathrm{s}}$.

$$\begin{array}{c}\hfill {\tilde{p}}_{\mathrm{s}}=\frac{2\pi }{x_{+}-x_{-}}\end{array}$$

(29)

The reason we set Equation (26) equal to 1, and the integral in Equation (22) equal to $2^N$, is because $\tilde{\mathrm{G}}$(x) is modeling the histogram of a QUBO’s solution space, like shown in Figure 2. This means that the total number of solutions to C(X) is $2^N$, and similarly the minimum number of distinct C(X${}_i$) solutions for a given cost function is 1. Therefore, after setting the integral in Equation (22) equal to $2^N$, solving $\tilde{\mathrm{G}}$(x)$=1$ yields approximations for C(X${}_{\min }$) and C(X${}_{\max }$) on the tails of the gaussian.

To demonstrate how well sampling is able to approximate Equation (14), we tested the strategy outlined above against the 1000 QUBOs from Figure 8 ($N=23$). For four values of M: 100, 500, 1000, and 2000, each QUBO was used for 50 trials of random sampling to produce approximate ${\tilde{p}}_{\mathrm{s}}$ values. These values were then compared to the true value of $p_{\mathrm{s}}$ from Equation (14), as given by Equation (30) below, and finally averaged together to produce Table 1.

$$\begin{array}{c}\hfill {\tilde{p}}_{\mathrm{s}}\mathrm{Error}=\frac{|{\tilde{p}}_{\mathrm{s}}-p_{\mathrm{s}}|}{p_{\mathrm{s}}}\end{array}$$

(30)

The significant result from Table 1 is that sampling 100 - 500 times, on a cost function of $2^{23}$ solutions, is accurate enough to produce an approximate ${\tilde{p}}_{\mathrm{s}}$ value with an expected error of only 7%. And as we show in the next section, this is enough accuracy to use for either a heuristic or variational approach for finding optimal solutions.

5. Variational Amplitude Amplification

The results of sections II - IV. demonstrate quantum’s aptitude for encoding and solving a QUBO problem using amplitude amplification. In this section we discuss how this potential can be realized from an experimental perspective. In particular, we focus on amplitude amplification’s ability to find optimal solutions under realistic circumstances with limited information. The results of this section are then used to motivate section VI., in which we discuss how amplitude amplification can be used in a hybrid classical-quantum model of computing, similar to other successful variational approaches [38,39,48].

5.1. Boosting Near-Optimal Solutions

The results shown in Figure 6, Figure 7 and Figure 8 focus on quantum’s potential for finding $|{\mathrm{X}}_{\min }\rangle$ and $|{\mathrm{X}}_{\max }\rangle$, the optimal solutions which minimize/maximize a given cost function C(X). However, in order to understand how amplitude amplification can be used in a variational model, it is equally as important that non-optimal $|{\mathrm{X}}_i\rangle$ states are also capable of boosting.

As discussed in the conclusion of our previous study [10], as well as sections III.C and IV.C, the most difficult aspect of using algorithm Algorithm 1 from an experimental standpoint is finding $p_{\mathrm{s}}$. More specifically, finding an optimal $p_{\mathrm{s}}$ for boosting $|{\mathrm{X}}_{\min }\rangle$ or $|{\mathrm{X}}_{\max }\rangle$ is a challenge due to the limited amount of information that one can learn through measurements alone. An example of this can be seen in Figure 9, which shows the peak achievable probabilities of the three lowest $|{\mathrm{X}}_i\rangle$ states as a function of $p_{\mathrm{s}}$ ($|{\mathrm{X}}_{\min }\rangle$ and the next two minimum solutions), for the QUBO corresponding to X${}_{\Delta }=331.5$ from Figure 7.

The challenge presented in Figure 9 is the narrow range of $p_{\mathrm{s}}$ values for which each $|{\mathrm{X}}_i\rangle$ state is able to achieve meaningful probabilities of measurement. From an experimental perspective, the $p_{\mathrm{s}}$ regions outside these bands are only capable of producing quantum superposition states which are slightly better than $|\mathrm{s}\rangle$, the equal superposition starting state. Thus, an experimenter could use a $p_{\mathrm{s}}$ value that is incredibly close to optimal, but only see seemingly random measurement results through repeat implementations of Algorithm 1.

Our proposed solution to the $p_{\mathrm{s}}$ problem as described above is twofold: 1) We must widen our view of useful $p_{\mathrm{s}}$ values and see where other $|{\mathrm{X}}_i\rangle$ states become highly probable, and 2) put less burden on quantum to find optimal solutions alone when an assisting classical approach may be more suitable. In this section we focus on addressing (1), which will then motive (2) in section VI.

Suppose we aren’t solely interested in using quantum to find the exact optimal solution C(X${}_{\min }$), but instead are content with any X${}_i$ within the best 50 answers (50 lowest C(X) values). In order for amplitude amplification to be viable in this heuristic context, it requires significant probabilities of measurement for these non-optimal solution states, similar to Figure 9. Additionally, an experimenter must be able to quickly and reliably find the $p_{\mathrm{s}}$ values which produce them. Shown below in Figure 10 is a plot which provides insight into the feasibility of both of these questions, for the QUBO corresponding to Figure 9.

Figure 10 shows the full $p_{\mathrm{s}}$ range for which an experimenter could find the 50 best solutions for minimizing C(X) via quantum measurements. The black circles indicate on the x-axis the $p_{\mathrm{s}}$ values where each $|{\mathrm{X}}_i\rangle$ state (or multiple states) is maximally probable, aligning with its corresponding C(X${}_i$) value along the y-axis. Numeric values for peak probabilities of the best 20 solutions are provided in the table below the plot, as well as a linear regression best fit (red-dotted line) for the overall 50 data points. The reported R correlation value is given by Equation (A5) in appendix $\mathbf{B}$.

There are several significant results displayed in Figure 10, the first of which requires returning to Equation (2). By limiting the allowed weighted values for $W_i$ and $w_{ij}$ to integers, all solutions to C(X) are consequently integers as well. This means that the linear correlation shown in the figure can in principle be used to predict $p_{\mathrm{s}}$ values where integer C(X${}_i$) solutions must exist. If $W_i$ and $w_{ij}$ are instead allowed to take on float values, which is more general of realistic optimization problems, the linearity of solutions like shown still persists but cannot be used for predictions of allowed C(X) values as reliably.

The linear best fit shown in Figure 10 is accurate for the top 50 solutions, but extending the $p_{\mathrm{s}}$ scale further reveals that it is only an approximation applicable to a small percentage of states. This is shown in Figure 11 below, which once again is a $p_{\mathrm{s}}$ vs. C(X) plot for the same QUBO, but now for the best 400 minimizing solutions. It is clear from the data in this figure that the top 400 solutions are in no way linearly aligned, which is a more expected result given the complicated nature of these imperfect gaussian distributions undergoing amplitude amplification. However, although the data is not linear, there is very clearly a curved structure that could be utilized to predict $p_{\mathrm{s}}$ values in the same manner described above.

It is important to note that in both Figure 10 and Figure 11, the manner in which the solution states $|{\mathrm{X}}_i\rangle$ are found to be most probable is sequential. This means that if a particular state $|{\mathrm{X}}_i\rangle$ is most probable at a certain value $p_{\mathrm{s}}=x$, all solutions C(X${}_j$) < C(X${}_i$) will have peak probabilities at values $p_{\mathrm{s}}<x$. However, the bottom plot in Figure 11 shows that the further a solution state is from $|{\mathrm{X}}_{\min }\rangle$ the lower its achievable peak probability. This means that there is a limit to how many solutions are viable for amplitude amplification to find. As we discuss in the coming subsections, these are the key underlying features that we must consider when constructing a variational amplitude amplification algorithm.

5.2. Constant Iterations

In order to construct an algorithm which capitalizes on the structure and probabilities shown in Figure 11, we must consider an additional piece of information not illustrated in the figure: step 3 of Algorithm 1, iterations k. The data points in the figure are indeed the $p_{\mathrm{s}}$ values which produce the highest probabilities of measurements, but unfortunately they are achieved using different iteration counts. In principle this means that an experimenter must decide both $p_{\mathrm{s}}$ and k before each amplitude amplification attempt, further complicating the information learned from measurement results.

Unlike $p_{\mathrm{s}}$, which is difficult to learn because it depends on the collective $2^N$ solutions to C(X), approximating a good iteration count k is easier. It turns out that the standard number of Grover iterations $k_{\mathrm{G}}=\frac{\pi }{4}\sqrt{N/M}$, where N is the total number of quantum states and M is the number of marked states, is equally applicable when using $U_{\mathrm{c}}$ as well. If an experimentor can use $k\approx k_{\mathrm{G}}$ iterations for a cost oracle $U_{\mathrm{c}}$ and find significant probabilities of measurment, then a variational amplitude amplification strategy can be reduced to a single parameter problem: $p_{\mathrm{s}}$. Figure 12 demonstrates that this is indeed viable, showcasing $|{\mathrm{X}}_i\rangle$ solution state probabilities as a function of $p_{\mathrm{s}}$ for three different choices of k.

The QUBO corresponding to Figure 12 is the same $N=25$ example for Figure 10 and Figure 11. For instances where multiple states correspond to the same numerical solution (C(X${}_i$) = C(X${}_j$)), the solid-color line shown represents their joint probability: Prob.( $|{\mathrm{X}}_i\rangle$ ) + Prob.( $|{\mathrm{X}}_j\rangle$ ) (note that these individual probabilities are always equal). Examples of this can also be seen in the table included in Figure 10. Additionally, a black-dashed line is shown in the top three plots, tracking the collective probability of the five most probable solutions at any given $p_{\mathrm{s}}$. These three lines are then replotted in the bottom panel for comparison.

The $p_{\mathrm{s}}$ region shown in Figure 12 was chosen to illustrate a scenario where variational amplitude amplification is most viable. For $p_{\mathrm{s}}>0.00291$, nearly every possible integer solution C(X${}_i$) $\ge -1497$ exists via some binary combination for this particular QUBO problem. The exceptions where certain integer solutions do not exist can be seen clearly in the $p_{\mathrm{s}}$ regions with very low probability, for example $0.0029065\le p_{\mathrm{s}}\le 0.002907$. Contrast to the region shown in this figure, once $p_{\mathrm{s}}$ becomes closer to where $|{\mathrm{X}}_{\min }\rangle$ is maximally probable, then measurment probabilities become more akin to Figure 9. Thus, it is more strategic for a hybrid algorithm to start in a $p_{\mathrm{s}}$ region like Figure 12, where measurement results can consistently yield useful information.

5.3. Information Through Measurements

From an experimental perspective, a significant result from Figure 12 are the black-dashed lines shown in the top three plots. At $k=3000$ ($k_{\mathrm{G}}\approx 4500$ for 25 qubits, $M=1$), the black-dashed line is almost entirely composed of the single most probable solution state(s). With probabilities around $70-80$% for many of the states shown, it is realistic that the same $|{\mathrm{X}}_i\rangle$ state could be measured twice in only $2-4$ amplitude amplification attempts. Two measurements yielding the same C(X${}_i$) value (possibly from different X${}_i$) is a strong experimental indicator that the $p_{\mathrm{s}}$ value used is very close to optimal for that solution, corresponding to the data points from Figure 10 and Figure 11. Confirming $3-4$ different data points in this manner can then be used to approximate the underlying curved structure of these figures, which in turn could be used to predict $p_{\mathrm{s}}$ values where $|{\mathrm{X}}_{\min }\rangle$ may exist.

While using k closer to $k_{\mathrm{G}}$ is good for getting the maximal probability out of solution states, the $k=1000$ and 2000 plots in Figure 12 support a different strategy for quantum. At $k=2000$, the black-dashed line is still primarily composed of the single most probable $|{\mathrm{X}}_i\rangle$ state(s), but critically it does not have the same dips in probability between neighboring solutions. Instead, the cumulative probability stays just as high for these in-between $p_{\mathrm{s}}$ regions, sometimes even higher! If we now look at the $k=1000$ plot, this trend becomes even more prevelant, whereby the cumulative probability plot is on average $20-30$% higher than any individual $|{\mathrm{X}}_i\rangle$ state. Interestingly, the bottom panel of Figure 12 shows that cumulative probability plot for $k=1000$ is higher than the $k=3000$ line in many regions. Thus, if the role of quantum is to simply provide a heuristic answer [49], not necessarily $|{\mathrm{X}}_{\min }\rangle$, then using lower k values is favorable for a few reasons. Firstly, we can anticipate solutions in a $p_{\mathrm{s}}$ region where multiple states share the same cost function value, so one can expect $M>1$ more frequently when using $k_{\mathrm{G}}=\frac{\pi }{4}\sqrt{N/M}$. Secondly, the amplitude amplification process itself is faster due to smaller k, which makes it more achievable on noisy qubits due to shallower circuit depths.

The optimal use of k is a non-trivial challenge to an experimentor. However, as illustrated in Figure 12, amplitude amplification can still be effective with a wide range of different k values. To further demonstrate this, Figure 13 shows three plots of simulated measurements over the $p_{\mathrm{s}}$ range depicted in Figure 12. Using the k values 1000, 2000, and 3000, each plot shows data points representing probabilistic measurements at regular intervals of $p_{\mathrm{s}}$. In order to compare the k value’s effectiveness more equally, the number of measurements taken per $p_{\mathrm{s}}$ value, t, was chosen such that $t·k=12000$ is consistent across all three experiments. Thus, each of the three plots in Figure 13 represents the same total number of amplitude amplification iterations divided among t experimental runs.

The data points shown in Figure 13 are separated into two categories, which are easily recognizable from an experimental perspective. Measurements which yielded C(X${}_i$) $<-1350$ are plotted as red circles, while all other measurements are plotted as black triangles. As illustrated for all three values of k, the red data points can be seen as producing near linear slopes, all of which would signal to the experimenter that these measurement results are leading to X${}_{\min }$. The motivation for Figure 13 is to demonstrate that the same underlying information can be experimentally realized using different k values. Thus, when to use $k=3000$ versus $k=1000$ is a matter of optimization, which we discuss in section VI. as the role of a classical optimizer for a hybrid model.

5.4. Quantum Verification

The results of the previous subsections demonstrate the capacity for amplitude amplification as a means for finding a range of optimal X${}_i$ solutions. However, regardless of whether these solutions are found via quantum or classical, a separate problem lies in verifying if a given solution is truly the global minimum X${}_i=$ X${}_{\min }$. If it is not, then X${}_i$ is refered to as a local minimum. Classically, evolutionary (or genetic) algorithms [50,51,52,53] are one example strategy for overcoming local minima. Similarly, quantum algorithms have also demonstrated success in this area for both annealing [54,55] and gate-based [56,57,58].

The strategy for verifying a local versus global minimum using amplitude amplification can be seen by comparing the region $0.0029\le p_{\mathrm{s}}\le 0.00291$ in Figure 12 and Figure 13. For the linear QUBO corresponding to these figures, there exists a solution C(X${}_i$) $=-1497$ which becomes maximally probable at $p_{\mathrm{s}}\approx 0.002914$, followed by the next lowest solution C(X${}_i$) $=-1491$ at $p_{\mathrm{s}}\approx 0.002892$. Because there are no binary combinations X${}_i$ that can produce values $-1492\ge$ C(X${}_i$) $\ge -1496$, the $p_{\mathrm{s}}$ region that $\mathit{would}$ correspond to their solutions instead produces nothing measurably significant. This can be seen by the low cumulative probabilities in Figure 12, as well as experimentally in Figure 13 as a gap in red data points for this $p_{\mathrm{s}}$ region across all three simulations.

The ability for quantum to determine if an X${}_i$ solution is locally or globally minimum is achieved by searching past the $p_{\mathrm{s}}$ value corresponding to the solution. Doing so will result in one of two outcomes: either a lower C(X${}_j$) value will be probabilistically found (confirming X${}_i$ was a local minimum), or the experimenter will only find random measurement results (X${}_i$ was the global minimum). Examples of this can be seen in Figure 14, showcasing simulated measurement results as an experimenter searches past the optimal $p_{\mathrm{s}}$ value for $|{\mathrm{X}}_{\min }\rangle$.

The simulated experiments shown in Figure 14 were chosen to highlight both favorable (bottom) and unfavorable (top) cases for quantum. The commonality between both experiments is that there is a clear point in $p_{\mathrm{s}}$ (grey line) in which decreasing $p_{\mathrm{s}}$ further results in only noisy random measurements. However, determining this cutoff point using measurement results alone is challenging. The top plot corresponds to the QUBO from Figure 10, Figure 11 and Figure 12, which is the non-ideal situation in which there are significant gaps in solutions between the best 20 minimizing C(X${}_i$). Experimentally this manifests as numerous $p_{\mathrm{s}}$ regions that could be wrongly interpreted as the X${}_{\min }$ cutoff point. Conversely, the bottom plot represents the ideal case where the best minimizing C(X${}_i$) solutions are all closely clustered together. This leads to a much more consistent correlation of measurement results leading to X${}_{\min }$, followed by an evident switch to randomness.

The significance here is that amplitude amplification has an experimentally verifiable means for identifying the global minimum X${}_{\min }$ of a cost function. Similarly, the same methodology can be in principle used to check for the existence of an X${}_i$ solution corresponding to any given cost function value, which we discuss further in section VII.C. However, the obvious drawback is that this verification technique relies on numerous amplitude amplification measurements finding nothing, which costs further runtime as well as being probabilistic. As we discuss in the next section, a more realistic application of this quantum feature is to help steer a classical algorithm past local minima, leaving the veification of X${}_{\min }$ as joint effort between both quantum and classical.

6. Hybrid Solving

The results of section V. were all features of amplitude amplification using $U_{\mathrm{c}}$ that were found through classical simulations of quantum systems. They represent the primary motivation of this study, which is to demonstrate amplitude amplifaction’s potential and the conditions for which it can be experimentally realized. By contrast, the discussions of section VI. here are more speculative. Given all of the results from sections III. - V., we now discuss how the strengths and weaknesses of amplitude amplification synergize with a parallel classical computer.

The plots shown in Figure 13 and Figure 14 represent a very non-optimal approach to finding X${}_{\min }$, functionally a quantum version of an exhaustive search. If the ultimate goal is to solve a cost function problem as quickly as possible, then it is in our best interest to use any and all tools available. This means using a quantum computer when it is advantageous, and similarly also recognizing when the use of a classical computer is more appropriate. In this section we discuss this interplay between quantum and classical, and the situations in which an experimenter may favor one or the other. Shown below in Figure 15 is the general outline of a variational amplitude amplification model which relies solely on quantum to produce X${}_{\min }$.

Given the current state of qubit technologies [59,60,61], performing one complete amplitude amplification circuit should be considered a scarce resource. As such, it is the role of a classical optimizer to determine the most effective use of this resource, choosing $p_{\mathrm{s}}$ and k values which will probabilistically get the most value out of each attempt. Determining optimal values to adjust a quantum circuit is the typical hybrid strategy found among other popular variational models of quantum computing [38,39,48]. The majority of the computational workload is placed on the QPU (quantum processing unit), while a classical optimizer is used in between runs to adjust quantum circuit parameters accordingly. As evidenced by Figure 13 and Figure 14, this model is possible for amplitude amplification as well. However, there is a different model of hybrid computing which better utilizes both quantum and classical’s strengths, shown below in Figure 16 and Figure 17.

The advantage to hybrid computing using the model shown in Figure 16 is that both processors are working in tandem to solve the same problem, utilizing information gained from one another. Information obtained through amplitude amplification measurements can be used to speedup a classical algorithm, and vice versa. As we discuss further in the next subsection, this pairing of quantum and classical is maximally advantageous when the strengths of both computers compliment each other’s weaknesses.

6.1. Supporting Greedy Algorithms

One notable strength of classical computing is `greedy’ algorithms, which oftentimes provide heuristic solutions for use cases ranging from biology and chemistry [49,62] to finance [63]. Greedy algorithms are particularly viable for problems that possess certain structures which can be exploited [64]. The key feature to these algorithms is that they focus on making locally optimal decisions which yield the maximal gain towards being optimal. Consequently, they are very good at finding near optimal solutions quickly, but are also prone to getting bottlenecked into local minima [65].

The motivation for pairing amplitude amplification with a classical greedy algorithm is best exemplified by Figure 12 and Figure 13. The quantum states illustrated in these figures represent $|{\mathrm{X}}_i\rangle$ states which rank as the ${30}^{\mathrm{th}}-{80}^{\mathrm{th}}$ best minimizing solutions to C(X). Under the right conditions it is reasonable to expect that a quantum computer could yield a solution in this range within $1-5$ amplitude amplification attempts. The question then becomes how quickly a classical greedy algorithm could achieve the same feat? Without problem specific structures to exploit, and as problem sizes scale like $2^N$, it becomes increasingly unlikely that classical can compete heuristically with quantum, which we argue is quantum’s first advantage over classical in a hybrid model.

Now, supposing that amplitude amplification does yield a low C(X${}_i$) solution faster than classical, the problem then flips back to being classically advantageous. This is because the X${}_i$ solution provided by quantum is now new information available to the classical greedy algorithm. As such, beginning the greedy approach from this new binary string is likely to yield even lower C(X${}_i$) solutions in a time frame faster than amplitude amplification. For example, this is the exact scenario in which genetic algorithms shine [50,51,52,53,63], where a near-optimal solution is provided from which they can manipulate and produce more solutions. And if a fast heuristic solution is all that is needed, then quantum’s job is done, and the best minimal solution found by the classical greedy algorithm completes the hybrid computation.

But if a heuristic solution is not enough, then we can continue to use a hybrid quantum-classical strategy for finding X${}_{\min }$. Referring back now to Figure 13 and Figure 14, the strategy for quantum is to use multiple amplitude amplification trials and measurements to approximate the underlying correlation from Figure 10 and Figure 11. The fastest means for achieving this is to work in a $p_{\mathrm{s}}$ region analogous to Figure 12, where experimentally one has the highest probabilities of measuring useful information. Simultaneously, the classical greedy algorithm can also find X${}_i$ solutions in this area as it searches for X${}_{\min }$. Knowledge of these X${}_i$ solutions can be directly fed back to quantum, which can be used to predict $p_{\mathrm{s}}$ values where solutions are known to exist, speeding up the process of determining a $p_{\mathrm{s}}$ vs. C(X) correlation. Thus, after quantum initially aided classical, subsequent information obtained from classical is then used to speed up quantum.

In the time it takes for quantum to experimentally verify enough $p_{\mathrm{s}}$ and C(X${}_i$) values to create a predictive correlation, we expect the classical algorithm to find a new lowest C(X${}_i$) solution, labeled X’${}_i$ in Figure 17. After investing additional trials into the amplitude amplification side of the computation, it is now time for quantum’s second advantage: verifying local versus global minima. Using an approximate $p_{\mathrm{s}}$ vs C(X) best-fit, the quantum computer can skip directly to the $p_{\mathrm{s}}$ value corresponding to best currently known X’${}_i$ solution. As discussed in section V.D, searching past this $p_{\mathrm{s}}$ value will result in one of two outcomes. Either the quantum computer will find a new best solution C(X${}_j$) < C(X’${}_i$), or confirm that X’${}_i$ is indeed the global minimum X${}_{\min }$. In the former case, the greedy algorithm now starts again from the new lowest solution X${}_j$, repeating this cycle between quantum and classical until X${}_{\min }$ is found. Figure 17 below shows a workflow outline of this hyrbid strategy.

The biggest advantage to using a hybrid model like shown in Figure 17 is that it can be adapted to each problem’s uniqueness. Problems with known fast heuristic techniques can lean on the classical side of the computation more heavily [66,67], while classically difficult problems can put more reliance on quantum [68,69]. But above all else, this model of computation incorporates and synergizes the best known classical algorithms with quantum, rather than competing against them.

7. More Oracle Problems

All of the results from sections III. - V. were derived from linear QUBOs according to Equations (1) - (4). However, these results can be applied to more challenging and realistic optimization problems provided that 1) all possible solutions can be encoded via phases by an appropriate oracle operation $U_{\mathrm{c}}$, and 2) the distribution of all possible answers is suitable for boosting the solution we seek (gaussian, polynomial, exponential, etc. [18]). Here we will briefly note some additional optimization problems which meet both of these criteria.

7.1. Weighted & Unweighted Max-Cut

The Maximum Cut problem (`Max-Cut’) is famously NP-Hard [68], where the objective is to partition every vertex in a graph $\mathbb{S}$ into two subsets ${\mathbb{P}}_1$ and ${\mathbb{P}}_2$ such that the number of edges between them is maximized. In the weighted Max-Cut version of the problem, each edge is given a weight $w_{ij}$, and the goal is to maximize the sum of weights contained on edges between ${\mathbb{P}}_1$ and ${\mathbb{P}}_2$. The unweighted Max-Cut problem has already been demonstrated as a viable use for amplitude amplification [17], which we will build upon further here via the weighted version. Equation (31) below is the cost function C(X) for the weighted Max-Cut problem, which can be converted to the unweighted case by setting every edge weight $w_{ij}=1$. The binary variables $x_i$ here represent being partitioned into ${\mathbb{P}}_1$ or ${\mathbb{P}}_2$.

Shown in Figure 18 is an example graph $\mathbb{S}$ and one of its solutions. This example graph is composed of 10 vertices, labeled 1 - 10, and a total of 15 connecting edges. Encoding this graph requires one qubit per vertex, where the basis states $|1\rangle$ and $|0\rangle$ represent belonging to the subsets ${\mathbb{P}}_1$ and ${\mathbb{P}}_2$ respectively. See the bottom graph in Figure 18 for an example solution state.

$$\begin{array}{ccc}\hfill \mathrm{C}\left(\mathrm{X}\right)& =& \sum _{\{i,j\}\in \mathbb{S}}{w}_{ij}|x_i-x_j|\hfill \end{array}$$

(31)

The cost oracle $U_{\mathrm{c}}$ for solving Max-Cut must correctly evaluate all $2^N$ solution states $|{\mathrm{X}}_i\rangle$ based on the edges of $\mathbb{S}$ according to Equation (31). For example, if vertices 1 and 2 are partitioned into different sets, then $U_{\mathrm{c}}$ needs to affect their combined states $|{\mathrm{Q}}_1{\mathrm{Q}}_2\rangle =|01\rangle$ and $|10\rangle$ with the correct phase, weighted or unweighted. Just like Figure 3 from earlier, we can achieve this with a control-phase gate CP($\theta$),with the intent of scaling by $p_{\mathrm{s}}$ later (see Figure 4). The caveat here is that we need this phase on $|01\rangle$ and $|10\rangle$, not $|11\rangle$, which means that additional X gates are required for the contruction of $U_{\mathrm{c}}$, shown below in Equation (32).

$$\begin{array}{ccc}\hfill \mathrm{X}& =& \left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]\hfill \end{array}$$

(32)

For the complete $U_{\mathrm{c}}$ quantum circuit which encodes the graph $\mathbb{S}$ in Figure 18, please see appendix $\mathbf{C}$. Once properly scaled by $p_{\mathrm{s}}$, the solutions which are capable of boosting are determined by the underlying solution space distribution of the problem, which can be seen in Figure 19 below. The histogram in this figure shows all $2^{10}$ C(X${}_i$) solutions to the graph $\mathbb{S}$ from Figure 18. Even for a 10 qubit problem size such as this, it is clear that the underlying solution space distribution shows gaussian-like structure.

One interesting feature of Max-Cut is that all solutions come in equal and opposite pairs. For example, the optimal solutions to $\mathbb{S}$ from Figure 19 are $|0100101110\rangle$ and $|1011010001\rangle$, which both yield 13 `cuts’. Mathematically there is no difference between swapping all vertices in ${\mathbb{P}}_1$ and ${\mathbb{P}}_2$, but physically it means that there are always two optimal solution states. Consequently, these states will always share the effect of amplitude amplification together, which is something an experimenter must be aware of when choosing iterations k.

Finally, moving from the unweighted to weighted version of Max-Cut increases the problem’s difficulty, but notably does not change the circuit depth of $U_{\mathrm{c}}$. Rather than using $\theta =1$ for all of the control-phase gates, each $\theta$ now corresponds to a weighted edge $w_{ij}$ of the graph. Similar to the QUBO distributions shown in Figure 7, this increase in complexity allows for more distinct C(X${}_i$) solutions, and consequently more variance in features such as ${\sigma }^{\prime }$ and X${}_{\Delta }$.

7.2. Graph Coloring

A similar optimization problem to Max-Cut is Graph Coloring, also known as Vertex Coloring [68], which extends the number of allowed partition sets ${\mathbb{P}}_i$ up to any integer number k ($k=2$ is equivalent to Max-Cut). Given a graph of vertices and edges $\mathbb{S}$, the goal is to assign every vertex to a set ${\mathbb{P}}_i$ such that the number of edges between vertices within the same sets is minimized. Shown below in Equation (33) is the cost function C(X) for a k-coloring problem, where the values of each vertex $x_i$ are no longer binary, but can take on k different integer values. The quantity inside the parentheses is equal to 1 if $x_i=x_j$, and 0 for all other combinations $x_i\ne x_j$. Just like with Max-Cut, setting all $w_{ij}=1$ is the unweighted version of the problem.

$$\begin{array}{ccc}\hfill \mathrm{C}\left(\mathrm{X}\right)& =& \sum _{\{i,j\}\in \mathbb{S}}{w}_{ij}\left(1-\lceil \frac{|x_i-x_j|}{k}\rceil \right)\hfill \end{array}$$

(33)

The name `coloring’ is in reference to the problem’s origins, whereby the sets ${\mathbb{P}}_i$ all represent different colors to be applied to a diagram, such as a map. Shown below in Figure 20 is an example picture composed of overlapping shapes, where each section must be assigned one of k colors such that the number of adjacent sections with the same color is minimized. Example solutions for $k=3$ and $k=4$ are shown, along with their vertex and quantum state representations of the problem.

In order to encode graph coloring as an oracle $U_{\mathrm{c}}$, the choice of k determines whether qubits or another form of quantum computational unit is appropriate. While qubits are capable of producing superposition between two quantum states, qudits are the generalized unit of quantum information capable of achieving superposition between d states [70,71,72,73]. To see why this is necessary, let us compare the $k=3$ and 4 examples from Figure 20, and the quantum states needed to represent partitioning each vertex.

For $k=4$, we need four distinct quantum states to represent a vertex belonging to one of the ${\mathbb{P}}_i$ partitions. While a single qubit can’t do this, a pair of qubits can. Thus, every vertex in $\mathbb{S}$ can be encoded as a pair of qubits, letting the basis states $|00\rangle$, $|01\rangle$, $|10\rangle$, and $|11\rangle$ each represent a different color. Alternatively, we could use a $d=4$ qudit to represent each vertex, assigning each partition a unique basis state: $|0\rangle$, $|1\rangle$, $|2\rangle$, or $|3\rangle$, such as the state shown in Figure 20. Mathematically the two encodings are identical, so the choice between whether to use qubits or qudits is a matter of experimental realization (i.e. which technology is easier to implement).

For $k=3$ however, two qubits is too many states, and a single qubit is not enough. So in order to represent three colors exactly in quantum, the appropriate unit is a `qutrit’ (the common name for a $d=3$ qudit). Similarly, all prime numbers d can only be encoded as their respective d-qudit, while all composite values can be built up from combinations of smaller qudits. Once an appropriate mixed-qudit quantum system is determined, constructing $U_{\mathrm{c}}$ is the same as the Max-Cut problem from earlier, but now with k state-state interactions. For an example of qudit quantum circuits and their use for amplitude amplification, please see Wang et al. [70] and our previous work on the Traveling Salesman problem[10].

7.3. Subset Sum

For all of the oracles discussed thus far, the circuit depth and total gate count for $U_{\mathrm{c}}$ is determined by the size and connection complexity of $\mathbb{S}$, the graphical representation of the problem. By contrast, the simplest possible quantum circuit that can be used as $U_{\mathrm{c}}$ corresponds to the Subset Sum problem [68]. The cost function for this problem is given in Equation (34).

$$\begin{array}{ccc}\hfill \mathrm{C}\left(\mathrm{X}\right)& =& \sum _i^N{W}_i{x}_i\hfill \end{array}$$

(34)

Rather than optimizing Equation (34), which is trivial, the Subset Sum problem is to determine if there exists a particular combination such that C(X${}_i$) = T, where T is some target sum value. The boolean variables $x_i$ represent which $W_i$ values to use as contributors to the sum. Figure 21 below shows an $N=10$ example. Note that this problem is equally applicable to any of the other oracles discussed thus far, whereby we can ask if a target value T exists for some graph $\mathbb{S}$.

The reason why Equation (34) is the simplest $U_{\mathrm{c}}$ oracle one can construct is because the cost function doesn’t contain any weights $w_{ij}$ that depend on two variables. Consequently, the construction of $U_{\mathrm{c}}$ doesn’t use any 2-qubit phase gates CP($\theta$), instead only requiring a single qubit phase gate P($\theta$) for every qubit. In principle, all of these single qubit operations can be applied in parallel, such as in Figure 3, which means that the circuit depth of $U_{\mathrm{c}}$ is exactly one.

Although this is the most gate efficient $U_{\mathrm{c}}$, using it to solve the Subset Sum problem comes with some limitations. Firstly, it can only solve for T values within a limited range. This is illustrated by the results of Figure 11, which demonstrate that amplitude amplification can only produce meaningful probabilities of measurement up to a certain threshold away from X${}_{\min }$ or X${}_{\max }$. Consequently, one can only use $U_{\mathrm{c}}$ here if the target sum value T is within this threshold distance from the extrema.

The second limitation to consider is the discussion from section V.D, whereby the information of whether a state C(X${}_i$) $=T$ exists or not may rely on measurements finding nothing. Previously we discussed how an experimenter might iteratively decrease $p_{\mathrm{s}}$ and eventually expect to find regions where cost function values do not exist (see Figure 14) as one approaches X${}_{\min }$. Here things are easier, since an experimenter can test for $p_{\mathrm{s}}$ values above and below where C(X${}_i$) $=T$ (except for the case where T is the global extrema). Using a $p_{\mathrm{s}}$ vs. C(X) correlation in this manner can confirm exactly where the $p_{\mathrm{s}}$ value for C(X${}_i$) $=T$ must be. Testing this $p_{\mathrm{s}}$ window will then either confirm the existence of a solution for T via a measurement, or conversely confirm no solution exists through multiple trials of random measurement results.

8. Conclusion

The results of this study demonstrate that the gate-based model of amplitude amplification is a viable means for solving combinatorial optimization problems, particularly QUBOs. The ability to encode information via phases and let the $2^N$ superposition of qubits naturally produce all possible combinations is a feature entirely unique to quantum. Harnessing this ability into a useful algorithmic form was the primary motivation for this study, and as we’ve shown, is not without its own set of challenges. In particular, the discussions of sections IV.A & IV.B highlight that this algorithm is not a `one size fits all’ strategy that can be blindly applied to any QUBO. Depending on how the numerical values of a given problem form a solution space distribution, it may simply be impossible for amplitude amplification to find one extrema or the other. Figure 8 shows that at least one of the extrema solutions is always viable for quantum to find, it just may not happen to be the one that is of interest to the experimenter.

For cases where the desired solution is well-suited for quantum to find, that is $|{\mathrm{X}}_{\min }\rangle$ or $|{\mathrm{X}}_{\max }\rangle$ is capable of achieving a high probability of measurement, a different challenge lies in finding the correct $p_{\mathrm{s}}$ value to use in order to boost these states. However, the results of section V. illustrate that this challenge is solvable via quantum measurement results. If the best an experimenter could do is simply guess at $p_{\mathrm{s}}$ and hope for success, then amplitude amplification would not be a practical algorithm. But the correlations shown in Figure 10 and Figure 11 illustrate that that is not the case, and that information about $p_{\mathrm{s}}$ can be experimentally learned and used to find extrema solutions. How quickly this information can be experimentally produced, analyzed, and used is exactly how quickly quantum can find the optimal solution, which is an open question for further research.

While the free parameter $p_{\mathrm{s}}$ can be considered the bottleneck of our algorithm for finding optimal solutions, there is a second important metric by which we can judge the usefulness of amplitude amplification: as a heuristic algorithm. A major finding of this study is depicted in Figure 12, which shows that there is a wide range of $p_{\mathrm{s}}$ values for which quantum can find an answer within the best $1-5$% of all solutions. And as we demonstrated in section IV.C with sampling, it is not unrealistic that a classical computation can estimate this $p_{\mathrm{s}}$ region very quickly. The question then becomes how does this compare to classical greedy algorithms, and how quickly can they achieve the same feat in a timescale compared to quantum’s O($\frac{\pi }{4}\sqrt{N/M}$) for problem sizes of $2^N$. The answer to this question will vary from problem to problem, but certainly in some cases such as highly interconnected QUBOs we view this as the first practical use for amplitude amplification.

And finally, there is one important sentiment from section VI that we would like to reiterate again here, namely that amplitude amplification is a technique that benefits tremendously from working in parallel with a classical computer. The information learned through quantum measurements can equally be of use to speeding up quantum as well as a classical algorithm. And vice versa, information learned through a classical greedy algorithm can be used to speed up quantum. The goal of this hyrbid computing model is to utilize the advantages both computers have to offer, and ultimately to find optimal solutions faster than either computer can achieve alone. Understanding which optimization problems this scenario may be applicable to is the future direction of our research.