How fast can quantum annealers count?

The paradigm of adiabatic quantum computation (AQC), proposed by Farhi et al [1] about a decade ago is a simple yet intriguing approach to problem solving on a quantum computer. Unlike the leading paradigm of circuit-based quantum computing, AQC is an analogue, time-continuous method which does not require the design and use of quantum gates. As such, it can in many ways be thought of as a simpler and perhaps more profound method for performing quantum computations while also being easier to design and implement experimentally [2–5].

The potential advantages of AQC over the usual circuit-based and other computational paradigms have generated substantial research, and intensive efforts have been made toward finding computational problems that are solved faster on an adiabatic quantum computer (henceforth, a quantum annealer) than on a classical computer (see, e.g., [6–15] and references therein). One of the most important results in that context, which establishes the universality and quantum power of the AQC paradigm is the existence of a theorem stating that adiabatic quantum computing is polynomially equivalent to circuit model quantum computation [16, 17].

Despite the above equivalence, it is not yet understood in which cases AQC is as efficient as circuit model quantum computation (or even potentially polynomially faster). It would therefore be interesting to look at problems for which circuit model algorithms are known, and to see whether equally-efficient adiabatic algorithms can be found. To date, however, there are almost no examples to show that AQC algorithms are similar in complexity to circuit-model algorithms. One notable exception is the adiabatic Grover's search algorithm devised by Roland and Cerf [18], demonstrating that the quadratic speedup over classical algorithms first shown by Grover [19] using circuit-based computation, can also be gained with adiabatic quantum algorithms (we discuss this algorithm in more detail later).

In what follows, we consider another problem for which, much like Grover's search problem, there is a circuit-model algorithm that is quadratically faster than its classical counterpart. This is the problem of 'quantum counting' [20, 21] in which one has to estimate the number of occurrences of a specific item in a database.

The problem of counting the number of solutions to an optimization problem (equivalently, the ground states of physical systems) is known to play a key role in the understanding of the physical properties of correlated classical as well as quantum systems (see, e.g., [22] and references therein). The computational difficulty of solving this problem is the focus of ongoing studies in computer science as well [22, 23].

Here, we show that the quantum counting problem, where one has to estimate the number of marked items in an unsorted 'blackbox' database, may also be solved on a quantum annealer and, quadratically faster than the best corresponding classical algorithms via the use of local adiabatic evolution (LAE) [18]. We shall see that the complexity of the proposed method is on par with the analogous circuit-based algorithm and so is quadratically faster than the corresponding classical algorithm.

The paper is organized as follows. In the next section we briefly discuss the principles of the quantum adiabatic algorithm (QAA) that is the heart of AQC. We then describe in some detail the concept of LAE in section 3. In section 4, we analyze the dynamical equations of the QAA for the adiabatic Grover's algorithm. These are later used in sections 5 and 6 to devise a QAA for solving the problem of quantum counting with quadratic speedup. Finally, in section 7, we conclude with a few comments.

We first describe the QAA which provides the general approach for solving problems by casting them into the form of optimization problems [1].

Within the framework of the QAA approach, the solution to an optimization problem is encoded in the ground state of a blackbox Hamiltonian $\hat{H}_p$. To find the solution, the QAA prescribes the following. As a first step, the system is prepared in the ground state of another Hamiltonian $\hat{H}_d$, commonly referred to as the driver Hamiltonian. The driver Hamiltonian is chosen such that it does not commute with the problem Hamiltonian and has a ground state that is fairly easy to prepare. As a next step, the Hamiltonian of the system is slowly modified from $\hat{H}_d$ to $\hat{H}_p$, using the linear interpolation, i.e.,

$$\begin{equation} \hat{H}(s)=s \hat{H}_p +(1-s) \hat{H}_d , \end{equation} \tag{ 1 }$$

where s(t) is a parameter varying smoothly with time, from 0 at t = 0 to 1 at the end of the algorithm, at $t=\mathcal {T}$. If this process is done slowly enough, the adiabatic theorem of quantum mechanics (see, e.g., [24] and [25]) ensures that the system will stay close to the ground state of the instantaneous Hamiltonian throughout the evolution, so that the system evolves to a state that is close to the ground state of $\hat{H}_p$. At this point, measuring the state will give the solution of the original problem with high probability.

In the simple case where s(t) varies from 0 to 1 at a constant rate, the runtime $\mathcal {T}$ must be chosen to be large enough so that the adiabatic approximation holds: this condition determines the efficiency, or complexity, of the QAA. A condition¹ on $\mathcal {T}$ can be given in terms of the instantaneous eigenstates {|m〉} and eigenvalues {E_m} of the Hamiltonian H(s), as [7, 30]

$$\begin{equation} \mathcal {T} >\frac{1}{\epsilon } \, {\max _{s} V_{01}(s) \over \min _s g^2(s)} , \end{equation} \tag{ 2 }$$

where is a small number, inversely proportional to the running time of the algorithm, that controls the expected probability for failure of the algorithm, g(s) is the first excitation gap E₁(s) − E₀(s) and V₀₁(s) = |〈0|dH/ds|1〉| (in our units ℏ = 1).

In the framework of LAE, the adiabatic parameter s(t) is varied not at a constant rate but rather at a variable rate, slowing down in the vicinity of the minimum gap and speeding-up in places where the gap is large [18]. The procedure is based on formulating a 'local' Landau–Zener condition for each value of the adiabatic parameter s. This condition can be written as:

$$\begin{eqnarray} &&\left| \frac{\mathrm{d}s }{\mathrm{d}t} \right| \le \epsilon \frac{g^2(s)}{V_{01}(s)} , \end{eqnarray} \tag{ 3 }$$

where g(s), V₀₁(s) and are as in (2).

When applicable, LAE can provide a significant speedup in calculation times [18, 31, 32]. A well-known example to demonstrate this, was first considered by Roland and Cerf [18]. In this example, an adiabatic algorithm for solving Grover's search problem [19] which describes a search for a marked item in an unstructured database, was found. More specifically, Roland and Cerf showed that while the application of the adiabatic theorem with a linear s(t) results in a running time that is of order N, where N is the number of items in the database, a carefully-chosen variable rate of s(t) yields a running time that scales with $\sqrt{N}$, i.e., a quadratic speedup is gained, similarly to the original result by Grover, that was found for the circuit-based model.

LAE however can only be efficiently used in cases where one has proper knowledge of the exact behavior of the gap and relevant matrix elements of the system for the problem at hand. This is normally not the case. The amenability of Grover's search problem to LAE stems from the fact that the gap of the system (as well as the relevant matrix element) can be calculated explicitly beforehand, even if the specific ground state is not known. With a proper choice of problem and driver Hamiltonians (which are discussed next), the gap in the Grover search problem is given by the simple expression:

$$\begin{eqnarray} &&g(\eta ;s)=[ 1-4(1-\eta )s(1-s) ]^{1/2} , \end{eqnarray} \tag{ 4 }$$

where η = M/N is the ratio of the number of solutions M to the size of the database N. Moreover, the matrix element $| \langle \mathrm{d}\hat{H} / \mathrm{d}s \rangle |$ can also be calculated in a straightforward manner to give:

$$\begin{eqnarray} &&V_{01}(s)=| \langle 0 | \mathrm{d}\hat{H} / \mathrm{d}s | 1 \rangle | = \frac{\sqrt{\eta (1-\eta )}}{g(\eta ; s)} . \end{eqnarray} \tag{ 5 }$$

The condition for LAE for the Grover problem is thus simply given by:

$$\begin{eqnarray} &&\left| \frac{\mathrm{d}s }{\mathrm{d}t} \right| \le \epsilon \frac{g^{3}(\eta ; s)}{\sqrt{\eta (1-\eta )}} . \end{eqnarray} \tag{ 6 }$$

A less strict inequality, using the fact that V₀₁(s) < 1, yields the more convenient bound of [18]:

$$\begin{eqnarray} &&\left| \frac{\mathrm{d}s }{\mathrm{d}t} \right| \le \epsilon g^{2}(\eta ; s) . \end{eqnarray} \tag{ 7 }$$

Integrating in time either of the two equations above between t = 0 (also, s = 0) and $t=\mathcal {T}$ (s = 1) yields in the limit of N → ∞:

$$\begin{eqnarray} &&\mathcal {T} \propto \frac{1}{\epsilon } \sqrt{N/M} = \frac{1}{\epsilon \sqrt{\eta }} \end{eqnarray} \tag{ 8 }$$

i.e., a quadratic speedup over the corresponding classical algorithms, a result also shared by Grover's original circuit-based algorithm [19].

An important property of the gap of Grover's search problem, (4), may be formulated as follows.

$$\begin{eqnarray} &&g(\eta _1;s) \le g(\eta _2;s) \quad \textrm {for any} \quad \eta _1 \le \eta _2 , \end{eqnarray} \tag{ 9 }$$

where η₁ and η₂ are two arbitrary ratios. The above property leads to the interesting (and somewhat surprising) attractive feature of LAE for the Grover search problem: if one runs the algorithm as if the ratio of solutions is η_* even though the actual ratio of solutions η may not be known, then one still obtains at the end of the adiabatic run a state that is very close the symmetric superposition of all solution states, as long as η_* ⩽ η. This is because choosing the adiabatic profile with η_*, namely,

$$\begin{eqnarray} &&\frac{\mathrm{d}s}{\mathrm{d}t} = \epsilon g^2(\eta _*; s) , \end{eqnarray} \tag{ 10 }$$

also satisfies the inequality (7) for any η_* ⩽ η. For example, a simple strategy for finding a solution state without knowing the number of solutions would be to run the algorithm with a profile function corresponding to M_* = 1 (i.e., η_* = 1/N). This will ensure a solution state at the end of the evolution with a (suboptimal) runtime of $\mathcal {T} \propto \frac{1}{\epsilon } \sqrt{N}$ (other, perhaps more efficient, strategies for finding a solution state without knowing M may also be found).

The above feature of the adiabatic procedure is to be contrasted with the circuit-based Grover's search algorithm, in which case one must know in advance the number of solutions in order to obtain a solution state with high probability. In the circuit-based case, the number of solutions M determines the number of 'Grover iterations' needed for the algorithm to end up in a solution state. Too many or too few iterations will drastically lower the probability of finding a solution. For that purpose, the quantum counting algorithm [20] has been devised, with the goal of estimating the number of solutions M prior to the application of Grover's search algorithm.

In the gate mode, the probability of success is periodic in the number of Grover iterations and it is this periodicity that has been utilized to construct the quantum counting algorithm. As discussed above, the probability of finding a solution in the adiabatic case is not periodic in the runtime of algorithm, so the techniques of the original quantum counting algorithm cannot, unfortunately, be used. However, in what follows we show that LAE, with the aid of the unique properties of the adiabatic Grover's problem presented above, may be used to solve the quantum counting problem adiabatically. The method suggested here, is comprised of repeated applications of the adiabatic algorithm for Grover's search problem of Roland and Cerf [18] but with varying, carefully-chosen, adiabatic profiles.

We therefore first analyze in what follows the dynamics of the adiabatic algorithm for Grover's search problem (a similar derivation is briefly discussed in [18]). Specifically, we derive here the equations describing an unsorted-database search with an unknown number (equivalently, ratio) of solutions M (or η) using an adiabatic profile s(t) with a ratio η_*. We then calculate the probability of a solution state at s = 1, in two limiting cases.

4.1. Dynamics of the adiabatic Grover's algorithm

Consider the problem Hamiltonian for the Grover problem [18]:

$$\begin{eqnarray} &&\hat{H}_p =1 - \sum _{m \in \mathcal {M}} | m\rangle \langle m | . \end{eqnarray} \tag{ 11 }$$

Here the M states $m \in \mathcal {M}$ are eigenstates of the computational basis with eigenvalue zero (all other states have eigenvalue one). The problem Hamiltonian is therefore a diagonal matrix of size N with M zero entries and (N − M) entries of one along the diagonal. The driver Hamiltonian is simply given by:

$$\begin{eqnarray} &&\hat{H}_d = - | \phi \rangle \langle \phi | , \end{eqnarray} \tag{ 12 }$$

where |ϕ〉 is the fully-symmetric state:

$$\begin{eqnarray} &&| \phi \rangle = \frac{1}{\sqrt{N}} \sum _{i=1}^N | i \rangle . \end{eqnarray} \tag{ 13 }$$

Obviously, the ground state of the driver Hamiltonian, i.e., the state that the system is initially prepared in, is

$$\begin{eqnarray} &&|\psi (t=0) \rangle = | \phi \rangle . \end{eqnarray} \tag{ 14 }$$

The Schrödinger equation

$$\begin{eqnarray} &&\mathrm{i}\frac{\mathrm{d}}{\mathrm{d}t} |\psi (t)\rangle =\hat{H} |\psi (t) \rangle \end{eqnarray} \tag{ 15 }$$

with |ψ(t)〉 = (c₀(t), c₁(t), ..., c_i(t)), ..., c_N(t)) simplifies considerably due to symmetry considerations: the state of the system has only two distinct amplitudes. The amplitudes c_i(t) that correspond to 'solution' states and those that correspond to 'non-solution' states. Let us denote the former set by a(t) and that of the latter by b(t). The Schrödinger equation thus becomes the two coupled equations:

$$\begin{eqnarray} &&\mathrm{i}\dot{a} = -(1-s) [\eta a+(1-\eta )b ] \end{eqnarray} \tag{ 16a }$$

$$\begin{eqnarray} &&\mathrm{i}\dot{b} = -(1-s) [\eta a+(1-\eta )b] +s b . \end{eqnarray} \tag{ 16b }$$

Note that N, the size of the problem and M, the number of solutions appear in the equations only as the ratio η = M/N.

The time-dependence of a and b is easily transformed into s dependence (as there is a one-to-one correspondence between the two in our case). The time-derivative transforms in the usual way: d/dt = (ds/dt) × d/ds, and the two coupled equations now become:

$$\begin{eqnarray} &&\mathrm{i}\epsilon g^2(\eta _*;s) a^{\prime } = -(1-s) [\eta a+(1-\eta )b] \end{eqnarray} \tag{ 17a }$$

$$\begin{eqnarray} &&\mathrm{i}\epsilon g^2(\eta _*;s) b^{\prime } = -(1-s) [\eta a+(1-\eta )b] +s b . \end{eqnarray} \tag{ 17b }$$

Here, the prime symbol (') stands for differentiation with respect to the new independent variable s and we have used (10) for ds/dt where the ratio η_* may in general be different than η. As already noted earlier, a choice of η_* ⩽ η results in a true adiabatic evolution and is expected to yield at the end of the QAA, a state that is very close to the equal superpositions of all solution states.

In addition, the adiabatic process dictates the initial conditions

$$\begin{eqnarray} &&a(s=0)=b(s=0)=1 , \end{eqnarray} \tag{ 18 }$$

which correspond to the fully-symmetric state (i.e., an equal superposition of all possible computational-basis states). The specific value of 1 chosen here leads to the simple normalization condition

$$\begin{eqnarray} &&\eta |a|^2+ (1-\eta ) |b|^2 = 1 , \end{eqnarray} \tag{ 19 }$$

that, similarly to the equations of motion, does not explicitly depend on N or M but only on their ratio.

4.2. End probabilities

At the end of the evolution, at s = 1, the final state of the system is a linear combination of the equal superposition of solutions states, namely, of $|\psi _{{\rm sol}} \rangle \propto \sum _{m \in \mathcal {M} } | m \rangle$ and the equal superposition of the non-solution states $|\psi _{{\rm no\hbox{-}sol}} \rangle \propto \sum _{m \notin \mathcal {M} } | m \rangle$. The closer the evolution is to being truly adiabatic, the larger the amplitude of the solution-states superposition will be. While the exact values of the two complex components at s = 1 cannot be computed analytically, one can obtain an expression for the probability of a solution state, i.e., the squared overlap between the final state at s = 1 and the optimal solution |ψ_sol〉:

$$\begin{eqnarray} &&P_{{\rm sol}} = |\langle \psi (s=1) | \psi _{{\rm sol}} \rangle |^2 = \lim _{s \rightarrow 1} \frac{\eta |a(s)|^2}{\eta |a(s)|^2+(1-\eta ) |b(s)|^2} , \end{eqnarray} \tag{ 20 }$$

in two important limiting cases. One is the limit of small , i.e., when the allowed error in the final state goes to zero. In the → 0 limit, the two equations are of the WKB type [33], and have a series solution of the form

$$\begin{eqnarray} &&a(s)=\mathrm{e}^{\frac{\mathrm{i}}{\epsilon } \sigma (s)} \quad \textrm {with}\quad \sigma (s)=\sum _{k=0}^{\infty } (\mathrm{i}\epsilon )^k \int ^s \alpha _k(s^{\prime })\ \mathrm{d}s^{\prime } , \end{eqnarray} \tag{ 21 }$$

where the functions α_k( · ) are real-valued. The above representation can be directly plugged into the coupled equations of motion. The first of the equations (17a) can be solved for b(s) to give:

$$\begin{eqnarray} &&b(s) = -\frac{\mathrm{e}^{\frac{\mathrm{i}}{\epsilon } \sigma (s)}}{1-\eta } \left(\eta - \frac{g^2 (\eta _*;s) }{1-s} \epsilon \sigma ^{\prime }(s) \right) \end{eqnarray} \tag{ 22 }$$

and the probability of a solution state simplifies to:

$$\begin{eqnarray} &&P_{{\rm sol}} = \lim _{s \rightarrow 1} \frac{\eta |a(s)|^2}{\eta |a(s)|^2+(1-\eta ) |b(s)|^2} = \lim _{s \rightarrow 1} \frac{\eta }{\eta +(1-\eta ) r(s)^2} , \end{eqnarray} \tag{ 23 }$$

where we have defined:

$$\begin{eqnarray} r(s)^2&=&|b(s)/a(s)|^2 \nonumber \\ &=& \left(\frac{\eta (1-s)-\epsilon g^2(\eta _*;s) \textrm { Re}[\sigma ^{\prime }(s)]}{(1-s)(1-\eta )} \right)^2 + \left( \frac{ \epsilon g^2(\eta _*;s) \textrm { Im}[\sigma ^{\prime }(s)]}{(1-s)(1-\eta )}\right)^2 . \end{eqnarray} \tag{ 24 }$$

Replacing σ(s) with its series expansion, turns the second equation (17b) into a series of equations for the various orders of . To second-order, we may keep only α₀( · ) and α₁( · ), at which stage solving for these and taking the s → 1 limit, yields:

$$\begin{eqnarray} &&P_{{\rm sol}}(\epsilon \rightarrow 0) = 1-\epsilon ^2 \eta (1-\eta ) +O(\epsilon ^4) . \end{eqnarray} \tag{ 25 }$$

A lengthier and more tedious calculation in which terms up to fourth-order of are kept, yields the fourth-order correction $\epsilon ^4 \eta (1-\eta ) ( 4-77 \eta +105 \eta ^2 +64 \eta _*-112 \eta \eta _* +16 \eta _*^2)$.

To leading order in , the probability to find a solution does not depend on value of η_*, as long as is small enough. As expected, when goes to 0 (i.e., when no error is allowed), the probability of the system to be in the equal superposition of solution states at the end of the QAA approaches 1 (in this case, the running time of the algorithm goes to infinity).

The more interesting limit however, is the one in which is small but kept fixed while N, the dimension of the Hilbert space (alternatively, the size of the database) goes to infinity. This limit translates to the limit where the ratios η = M/N and η_* = M_*/N are small², i.e., M, M_* ≪ N. In this case, the probability of success is completely determined by the Landau–Zener transition in the immediate vicinity of the minimum gap, at s = 1/2. There, the gap can easily be calculated to be $g(\eta ;s=1/2)=\sqrt{\eta }$. In addition, the Landau–Zener sweep rate is determined on the other hand by dg/dt = dg/ds × ds/dt = 2 × g²(η_*, s = 1/2) = 2η_*. The Landau–Zener formula therefore immediately yields:

$$\begin{eqnarray} &&P_{{\rm sol}}(N \rightarrow \infty ) \approx 1-\mathrm{e}^{-\frac{\pi \eta }{4 \epsilon \eta _*}} = 1 - \mathrm{e}^{-\frac{\pi M}{4 \epsilon M_* }} . \end{eqnarray} \tag{ 26 }$$

As expected, the probability of obtaining a solution increases with decreasing and also, as discussed in the previous section, with decreasing M_*. Most importantly, in the N → ∞ limit, the above probability has a well established N-independent value. The above probability will become useful in what follows when we address the problem of quantum counting.

Before addressing the full-complexity of the general quantum counting problem, let us first consider the special case of finding whether or not the number of solutions M (number of zero entries in the problem Hamiltonian) is zero. Circuit-based algorithms give a running time that is on the order of $O(\sqrt{N})$ [20, 21], that is, a quadratic speedup over the analogous classical algorithm.

To solve this problem quantum-adiabatically, one executes the QAA algorithm with the adiabatic profile function s(t) that obeys ds/dt = g²(η = 1/N, s), i.e., a profile that is suitable for solving the case of any nonzero number of solutions with quantum speedup.

Now, if the number of solutions M is zero, then the problem Hamiltonian will have no effect on the state of the system, and the state will remain the initial state |ϕ〉 throughout the adiabatic evolution, except for an undetectable global phase. A measurement of the z-magnetization at the end of the run will produce a random configuration with energy 1.

If on the other hand the number of solutions is nonzero, the Grover LAE procedure ensures us that at the end of each of the run, at s = 1, the final state of the system will be very close to |ψ_sol〉, the equal superposition of zero-energy states with a probability greater than $P_{{\rm sol}} = 1 - \mathrm{e}^{-\frac{\pi M}{4 \epsilon }}$ (in the N → ∞ limit), and a measurement along the z-direction will produce a zero-energy state with this probability.

We conclude then, that one execution of the adiabatic algorithm with runtime $\mathcal {T} \propto \sqrt{N}/\epsilon$ is sufficient for establishing whether or not the number of solutions is zero.

The general case of estimating the number of solutions M can be analyzed with the aid of the probability P_sol obtained in (26). Let us first briefly analyze the classical analogue of that same problem. In the classical case, random sampling of items in the database will give $P_{{\rm sol}}^{{\rm (classical)}}=M/N=\eta$ for each picked item. The number of successes after k trials will therefore be distributed binomially with the unbiased mean of $\hat{\eta }^{{\rm (classical)}} = \eta$ and a sampling error of:

$$\begin{eqnarray} &&\Delta \hat{\eta }^{{\rm (classical)}} \sim \sqrt{\frac{\eta (1-\eta )}{k}} \approx \sqrt{\frac{\eta }{k}} . \end{eqnarray} \tag{ 27 }$$

In terms of the number of solutions M, the error in M is thus

$$\begin{eqnarray} &&\Delta M^{{\rm (classical)}} \sim \sqrt{\frac{M N}{k}} , \end{eqnarray} \tag{ 28 }$$

i.e., the number of trials k required for estimating the number of solutions is

$$\begin{eqnarray} &&k \sim \frac{N M}{\Delta M^2} . \end{eqnarray} \tag{ 29 }$$

Specifically, for an error estimate of $\Delta M \sim \sqrt{M}$, the number of trials k required, has to scale linearly with the dimension of the database N.

In the quantum-adiabatic case, the question thus becomes: how many runs k of the adiabatic Grover's search algorithm (and with which adiabatic profiles), each run producing a solution state with probability P_sol [as given by (26)], are necessary to yield an unbiased estimate of M?

We answer this question by proposing the following two-step algorithm: first, we find a suitable value of M_*, according to which the adiabatic paths for subsequent QAA runs will be constructed. As a second step, we calculate how many QAA runs, with that chosen value of M_* are required to estimate M for a given error rate ³.

The first step of the algorithm is designed to find an optimal value for M_* that will determine the adiabatic profile to be used in the second step. Note that M_* should be on the order of the unknown M. Choosing a value for M_* that is too large (i.e., M_* ≫ M/) would correspond to probabilities that are very close to zero, leading to very poor statistics that would in turn not suffice for the determination of M in an accurate manner. Choosing a value for M_* that is too small (i.e., M_* ≪ M/) would lead to unnecessarily-long running times for the adiabatic runs and would correspond to probabilities that are very close to one, leading again to very poor statistics. An optimal choice for M_* would thus be a value that would yield a significantly-large probability for obtaining a solution state (one reasonable choice could be, e.g., 1 − P_sol ≈ 0.1), which implies M_* ∼ M/.

Since M is unknown, a way to find a value of M_* that would yield a given probability, is to employ a simple noisy binary search algorithm [34–36]. Starting with geometrically diminishing values of M_*, namely, M_* = N/2, N/4, ..., the QAA should be executed multiple times for each value, such that for each tested value M_*, the average ratio of solutions is computed (note that since M ≪ N the ratios will be zero at first). When M_* comes close to the vicinity of M and the obtained ratios of solutions are no longer strictly zero, the trial M_* values will need to be adjusted (i.e., increased or decreased) according to rules of binary search, until an approximate value of the chosen probability of solutions, P_sol is obtained. Note that in order to obtain a certain probability of solutions P_sol with a square-root error accuracy, only a fixed (i.e., an M-independent) number of QAA executions is required for each tested value of M_* (as the number of 'solution' outcomes is distributed binomially). For a binary search such as the one outlined above, the number of trial M_* values required is logarithmic in N, where the complexity of each trial run is $O(\sqrt{N/M_*})$ [as is given in (8)]. The overall complexity of this first step is therefore simply $O(\sqrt{N/M})$.

Having found an optimal value for M_*, the second step of the algorithm is initiated, consisting of the execution of multiple QAA runs with the chosen value of M_* (note that the chosen value obeys M_* ∼ M/). A straightforward calculation analogous to the one performed for the classical case above, following (27), yields, in the quantum-adiabatic case, the number of trials:

$$\begin{eqnarray} &&k \propto \frac{16 \epsilon ^2 M_*^2 \Big(\mathrm{e}^{\frac{\pi M}{4 \epsilon M_*}}-1\Big)}{\pi ^2 (\Delta M)^2} \approx \frac{4 \epsilon M_* M}{\pi (\Delta M)^2} , \end{eqnarray} \tag{ 30 }$$

where the last approximation is justified by the smallness of the probability 1 − P_sol. Since the complexity of each QAA run is $O(1/\epsilon \sqrt{N/M_*})$, the overall complexity of the algorithm will be $O(k/\epsilon \sqrt{N/M_*}) = O((M/\Delta M)^2 \sqrt{ N/\epsilon M})$. For an error that scales with the square root of the number of solutions, i.e., $\Delta M \sim \sqrt{M}$, the complexity of the algorithm will be $O(\sqrt{N M/\epsilon })$, i.e., the quadratic speedup over the classical analogue is obtained, as the complexity only scales with the square root of the size of the database. Comparing to the corresponding quantum circuit algorithm [20], we find that while both algorithms scale similarly with the size of the database, the scaling of M is somewhat different. For example, for an $O(\sqrt{M})$ error, the circuit-based algorithm yields an $O(\sqrt{N})$ complexity, i.e., it is $\sqrt{M}$ faster than the QAA proposed here. Since, however, we are considering here the case of small number of solutions, this difference plays only a minor role.

A similar comparison for the case where the error is required to scale linearly with the number of solutions, i.e., ΔM ∼ M, the adiabatic algorithm above yields a complexity of $O(\sqrt{N/\epsilon M})$, which in this case is identical to the circuit-based scaling [20]. For the convenience of the reader, we have added a table comparing the complexities of the classical, quantum-adiabatic and quantum circuit algorithms for different error estimates (see table 1).

Table 1. Complexities of the classical, quantum circuit and quantum-adiabatic counting algorithms for different error estimations. Here, N is the size of the database, M is the number of solutions and ΔM is the error estimation.

Error scaling	Classical algorithm	Circuit quantum algorithm [20]	Quantum-adiabatic algorithm
ΔM ∼ O(1)	O(NM)	$O(\sqrt{N M})$	$O(M \sqrt{N M})$
$\Delta M \sim O(\sqrt{M})$	O(N)	$O(\sqrt{N})$	$O(\sqrt{N M})$
ΔM ∼ O(M)	O(N/M)	$O(\sqrt{N/M})$	$O(\sqrt{N/M})$

We have given a prescription for solving the quantum counting problem in which the number of occurrences of a specific item within an unsorted database is to be estimated. Using local adiabatic evolution on Grover-type problems, the algorithm was shown to perform quadratically faster than the corresponding classical algorithm and on par with the circuit-based result. We have thus shown that quantum-adiabatic evolution may be utilized toward counting the number of solutions (or ground states) to a given problem and is not restricted to only finding those.

We note here, that similarly to the original Grover's algorithm, the problem Hamiltonian of the quantum counting algorithm presented here, is given as an oracle, and the quadratic speedup exhibited here holds only relative to it.

It is known that the gate-model quantum counting algorithm, may be viewed as a special case of a quantum phase estimation algorithm [21], in which one is required to estimate an eigenvalue of an implementable unitary operator. It would be interesting to see whether the algorithm outlined above can be used as the basis for a similar adiabatic calculation. Since however in the adiabatic case one deals with Hermitian operators, as opposed to the unitary operators of the circuit model, it is unclear how this could be done. Nonetheless, it is worth noting that the fully-symmetric initial wave function |ϕ〉 may be re-expressed in the form [21]:

$$\begin{eqnarray} &&| \phi \rangle = \sqrt{\eta } | \alpha \rangle + \sqrt{1-\eta } | \beta \rangle , \end{eqnarray} \tag{ 31 }$$

where we have defined:

$$\begin{eqnarray} &&| \alpha \rangle = \frac{1}{\sqrt{N \eta }} \sum _{m \in \mathcal {M}} | m \rangle , \end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray} &&| \beta \rangle = \frac{1}{\sqrt{N(1-\eta )}} \sum _{m \notin \mathcal {M}} | m \rangle . \end{eqnarray} \tag{ 33 }$$

Replacing η = M/N with the appropriate angle of sin ²(θ/2), we can rewrite the wave function as:

$$\begin{eqnarray} &&| \phi \rangle =\sin \left( \theta /2\right) | \alpha \rangle + \cos \left( \theta /2\right) | \beta \rangle , \end{eqnarray} \tag{ 34 }$$

which reveals that estimating M is equivalent to estimating an 'angle' θ that is associated with the problem Hamiltonian, which could possibly be linked to a phase of a unitary operator.

The importance of the above results, aside from the implications resulting from having a quadratically-faster quantum-adiabatic algorithm for this specific problem, lies in the fact that it is one of very few examples that illustrate the power and potential encompassed in adiabatic quantum computation, providing an adiabatic quantum algorithm that is as efficient as its circuit-model counterpart.

We thank Peter Young, Eleanor Rieffel and Jérémie Roland for useful comments and discussions. We acknowledge partial support by the National Security Agency (NSA) under Army Research Office (ARO) contract number W911NF-09-1-0391, and in part by the National Science Foundation under grant no. DMR-0906366.