Last time we looked at the basic theory behind quantum search based on the Grover’s algorithm. We went through the most basic case, a data set consisting of four items, and applied the algorithm to that, learning in the process that it managed to find the relevant entry we were looking for in a single step - compared to an average expected 2.25 steps required by the classical computation theory.

In this part, we will look at the more general theory behind Grover’s algorithm, and implement the general Q# variant that can be used to find any number in an arbitrarily large data set.

Oracles, again 🔗

Just like this was the case in our earlier discussions on the basic Deutsch’s problem and Deutsch-Jozsa algorithms, the theoretical model for a general Grover’s algorithm is based around the concept of the oracle function, used to determine whether the solution to the problem. In fact there is an entire class of such oracle based problems in quantum information theory.

In our case - the problem of finding a proverbial needle in haystack, or rather a single result in an unsorted data set - the oracle is a verification function $f(x)$, which determines whether we found the expected result or not. Suppose we want to find $y$. The verification function then returns 1 (true) or 0 (false) in the following way:

$$
f(x) =
\begin{cases}
1 & \text{when } x = y \
0 & \text{when } x \neq y \
\end{cases}
$$

Additionally, the oracle will need to use an extra ancilla qubit called the “oracle qubit”. The presence of the oracle qubit allows reversibility of the computation, which is mandated by the laws of quantum mechanics. Because of the presence of the ancilla qubit, the oracle’s transformation can be summarized as:

$$
\ket{x}\ket{y} = \ket{x}\ket{y \oplus f(x)}
$$

where

$$\ket{x} = \ket{0^{\otimes n}}$$

In other words, the input into the oracle can be an arbitrary $n$ qubit state $\ket{x}$, and the single oracle qubit $\ket{y}$ will get flipped if $f(x)$ returned true.

Grover’s algorithm theory 🔗

The construction of Grover’s algorithm resembles, or indeed matches, in its initial phases, Deutsch-Jozsa algorithm. As such, the preparation step of the algorithm will include creating a uniform superposition of all the input qubits. Additionally, we will rely here on the exact same phase kickback trick that we used in the Deutsch-Jozsa algorithm, allowing us to convert the change of sign on the oracle qubit to a phase shift of $\ket{x}$. For the phase kickback to work the oracle qubit needs to be in the initial state $\ket{1}$ instead of the usual $\ket{0}$, and it needs to be in the superposition as well.

We can describe the entire preparatory procedure in the following algebraic transformation chain, where $\ket{x}$ is an abritrary number of qubits and $\ket{1}$ is the ancilla qubit:

$$
\ket{\Psi_1} = (W \otimes H)\ket{x}\ket{1}
$$

As we learnt already in the earlier parts, this creates the following state:

$$\ket{\Psi_1} = \frac{1}{\sqrt{2^n}}\sum_{x \in {0,1}^n}\ket{x}\ket{-}
$$

It is common to use $2^n = N$, and because of that, we can rewrite $\ket{\Psi_1}$ into:

$$\ket{\Psi_1} = \frac{1}{\sqrt{N}}\sum_{x=0}^{N-1}\ket{x}\ket{-}
$$

This provides us with the foundational pieces for us to work with - the input into the oracle, which is called next. If use the notation $U_f$ for the oracle, keeping in mind that it realizes the $f(x)$ function define above, we arrive at state $\ket{\Psi_2}$:

$$\ket{\Psi_2} = U_f\ket{\Psi_1}
$$

This leads us to the following state, including the phase kickback, which we already know from the Deutsch-Jozsa algorithm:

$$\ket{\Psi_2} = \frac{1}{\sqrt{N}}\sum_{x = 0}^{N-1}(-1)^{f(x)}\ket{x}\ket{-}
$$

Thanks to the phase kickback, which is really the consequence of the linearity of mathematical formalism employed by quantum mechanics, the global sign changes to a minus. The oracle is therefore said to “mark” the solution for the search problem, as it changes the global phase of our quantum state. A side effect of that is that the oracle qubit actually remains unchanged, since its flip ends up being converted into the global phase shift, and stays in the state $\ket{-}$. As such, we will actually be able to ignore it from this point on, as it is unentangled from the rest of the computation.

While the marking of the answer using the phase is very helpful, as we learnt in the last part, the probability amplitudes are still unchanged, meaning we are not closer to reading out the correct solution yet. Therefore, the next step is to perform the inversion about the mean, which will allow us to amplify the probability amplitude of the marked solution.

To do that, let’s consider how can we formalize the inversion about the mean procedure. We discussed it last time, with a simple example of four probability amplitudes - three of them equal to $\frac{1}{4}$ and one with a phase shift applied, having the opposite sign $-\frac{1}{4}$. Let’s generalize this mathematically now. If we consider $m$ to be the mean in some larger data set and $v$ to be some value we will want to be inverted about the mean, the formula to obtain the inverted result $v’$ can be expressed as:

$$v’ = m + (m - v) = 2m - v
$$

We have to emphasize that in the quantum version of our search problem, we deal not with a sequence of loose numbers, but rather a single multi dimensional vector. Therefore, we can rephrase the above formula into the following algebraic representation, where $M$ is the matrix that finds the mean of the vector’s components, $I$ is the identity matrix and $V$ and $V’$ are multi dimensional vectors.

$$V’ = V(2M - I)
$$

We are now in a position to define the next step of the algorithm, the so-called Grover diffusion operator. It can be realized as the conditional phase shift $P$, surrounded by $W$ ($H^{\otimes n}$) transformations. The conditional phase shift will leave the sign of $\ket{x}$ as-is if it equals to $0^n$ and flip the sign in all other cases.

$$
P = \ket{x} \rightarrow
\begin{cases}
\ket{x} & \text{when } x = 0^n\
-\ket{x} & \text{otherwise}
\end{cases}
$$

We can therefore write the definition of $P$ as:

$$
P = 2\ket{0^n}\bra{0^n} - I
$$

The $\ket{a}\bra{b}$ syntax used here is the algebraic outer product. By applying a $W$ transformation before and after $P$ we obtain the final definition of the diffusion operator $D$:

$$WPW = 2\ket{\varphi}\bra{\varphi} - I = D
$$

where $\varphi$ denotes a uniform superposition of all the input states. For our theoretical model of Grover’s algorithm $\ket{\varphi}$ is represented by the state $\ket{\Psi_2}$ - the state after the application of the oracle $U_f$. We already are in a uniform superposition at that point, with one of the items having a phase changed. This allows us to formulate state $\ket{\Psi_3}$:

$$\ket{\Psi_3} = D\ket{\Psi_2}
$$

The diffusion operator, as already stated, inverts the probabilities around the mean, thus amplifying them, and increasing our classical probabilities of reading out a correct result.

Together, $D$ and $U_f$, make up the so called Grover iteration, $G$.

$$G = (2\ket{\varphi}\bra{\varphi} - I)U_f = DU_f
$$

In the simple 4 element data set example from last post, a single Grover iteration was already enough to obtain the answer with probability equal to 1 - we could proceed to the measurement process. However, in case of larger data sets, we need to repeat the Grover iteration about $\sqrt{N}$ times. Subsequent iteration, and applications of both the oracle $U_f$ and the diffusion $D$, amplify the probabilities further - up to a certain point, after which the probability starts to degrade and we cannot amplify further beyond that. To be more precise, that process is cyclical and after certain amount of amplifications it will start improving the results again. There is a geometrical way of calculating precisely the most appropriate amount of Grover iterations, though studying this is beyond the scope of this post.

To summarize, after $G$ is applied sufficient amount of times, we can perform the measurement on $\ket{x}$, which should give us the correct result with a very high degree of probability. Let’s recap the entire procedure:

1. Prepare state $\ket{00..0}$
2. Apply W transformation to create a uniform superposition
3. Apply the marking operation (oracle) $U_f$ using the extra ancilla qubit prepared in state $H\ket{1}$
4. Apply the inversion about mean (diffusion) $D$
5. Repeat steps 3 and 4 about $\sqrt{N}$ times
6. Measure the qubits to obtain a result with very high probability

Q# implementation 🔗

In this section, equipped with the newly acquired theoretical knowledge from the theoretical part of this post, we will go through a more advanced Grover’s search implementation than what we did in the last article. However, we shall use some part of that code here.

We shall start by introducing a general outline for a Q# operation that will perform the Grover search across an arbitrary number of qubits, for an arbitrary marked result. The snippet below shows such high level structure which prepares a uniform superposition of the initial state, and contains the necessary measurement already using the $LittleEndian$ register, and a direct conversion into an integer using $MeasureInteger$ operation from the $Microsoft.Quantum.Arithmetic$ namespace. We also define the amount of iterations as something that can be manually defined, because, as we said, an accurate calculation of this is not in scope of this post.

operation Grover(markIndex : Int, numberOfQubits : Int, iterations : Int) : Int {
    use qubits = Qubit[numberOfQubits];
    
    // apply W
    ApplyToEachA(H, qubits);
    
    // implement Grover iterations
    
    // measure
    let register = LittleEndian(qubits);
    let number = MeasureInteger(register);
    ResetAll(qubits);

    return number;
}

What is left is the implementation of the Grover iteration - both the “marking” part that will act as an oracle, and the amplitude amplification. For now, for simplicity, we can implement all of that inline in a single loop definition.

In the last post, we used the $CZ$ transformation to flip the sign (“mark”) the phase of $\ket{11}$ state, and then we manually moved the phase onto the number that interested us by applying the relevant $X$ transformations. This was easy enough to be done using a set of if statements, because we only dealt with 2 bit integers. Here, since we support an arbitrary number of bits (and thus, qubits), we need to generalize this solution. To do that we can implement the following procedure:

1. Take the integer that we want to mark (“the solution”) and produce a binary representation of it - a boolean array, that has the same length as the length of our qubit register. Let’s call this “solution binary array”
2. Use the solution binary array to apply $X$ transformation to our quantum register at every index in which the solution binary array has a bit value 0. This creates for us the state $\ket{11..1}$
3. Apply $CZ$ to flip the phase of it
4. Apply the same computation as in step 2. This will undo that operation and return us to the state before $CZ$ was applied, with the phase flip of $CZ$ moved onto the solution

After that we can apply the amplitude amplification as we did in the previous post, and we can repeat the entire process the required number of times - namely $\sqrt{N}$ times.

The entire Grover iteration written linearly using the above recipe is shown in the snippet below:

for x in 1..iterations {

    let markerBits = IntAsBoolArray(markIndex, numberOfQubits);
    for i in 0..numberOfQubits-1
    {
        if not markerBits[i] {
            X(qubits[i]);
        }
    }
    
    Controlled Z(Most(qubits), Tail(qubits));

    for i in 0..numberOfQubits-1
    {
        if not markerBits[i] {
            X(qubits[i]);
        }
    }

    // amplitude amplification
    ApplyToEachA(H, qubits);
    ApplyToEachA(X, qubits);
    Controlled Z(Most(qubits), Tail(qubits));
    ApplyToEachA(X, qubits);
    ApplyToEachA(H, qubits);
}

This snippet can now be inserted into the earlier code example, to produce the final version of the operation:

operation Grover(markIndex : Int, numberOfQubits : Int) : Int {
    use qubits = Qubit[numberOfQubits];

    // superposition
    ApplyToEachA(H, qubits);

     // Grover iterations
    for x in 1..iterations {

         // mark the solution
        let markerBits = IntAsBoolArray(markIndex, numberOfQubits);
        for i in 0..numberOfQubits-1
        {
            if not markerBits[i] {
                X(qubits[i]);
            }
        }
        
        Controlled Z(Most(qubits), Tail(qubits));

        for i in 0..numberOfQubits-1
        {
            if not markerBits[i] {
                X(qubits[i]);
            }
        }

        // amplitude amplification
        ApplyToEachA(H, qubits);
        ApplyToEachA(X, qubits);
        Controlled Z(Most(qubits), Tail(qubits));
        ApplyToEachA(X, qubits);
        ApplyToEachA(H, qubits);
    }

    let register = LittleEndian(qubits);
    let number = MeasureInteger(register);
    ResetAll(qubits);

    return number;
}

We can test whether the solution works, by implementing the entry point to our Q# program that invokes the operation and passes in the relevant arguments to configure the Grover search. For example, the following snippet runs the search on four qubits - where the data set ranges from 0 to 15, with 10 as the marked number.

@EntryPoint()
operation Main() : Unit {
    let result = Grover(10, 4, 3);
    Message($"Measured Grover example: {result}");
}

While we won’t calculate it here, after 3rd iteration, the probability of reading out 10 equals 96.1319% and we are very likely to read out 10 and see:

Measured Grover example: 10

If we increase the data set to 10 qubits, the data range grows to 1024 values.

@EntryPoint()
operation Main() : Unit {
    let result = Grover(10, 10, 25);
    Message($"Measured Grover example: {result}");
}

In that case $\sqrt{N} = 32$, but the most optimal amount of iterations is actually 25, when the probability of reading our marked value is 99.9461%.

Measured Grover example: 10

Finally, we can improve our code a bit by taking advantage of some of the features of Q#. First of all, both the marking code and the amplitude amplification could be extracted to standalone operations to reduce cluttering.

First let’s do that with the amplitude amplification code. Here we can additionally take advantage of the conjugation feature of Q#. If we have the transformation $V$ which is surrounded by transformations $U^\dagger$ and $U$, or rather a specific transformation $U$ and its reverse, we can use the $apply…within$ syntax of Q# for our convenience. Thus, instead of writing:

ApplyToEachA(H, qubits);
ApplyToEachA(X, qubits);
Controlled Z(Most(qubits), Tail(qubits));
ApplyToEachA(X, qubits);
ApplyToEachA(H, qubits);

We can write:

within {
    ApplyToEachA(H, qubits);
    ApplyToEachA(X, qubits);
} apply {
    Controlled Z(Most(qubits), Tail(qubits));
}

Therefore, our amplification operation will look like this:

operation AmplifyAmplitude(qubits : Qubit[]) : Unit is Adj {
    within {
        ApplyToEachA(H, qubits);
        ApplyToEachA(X, qubits);
    } apply {
        Controlled Z(Most(qubits), Tail(qubits));
    }
}

Additionally, turns out that Q# already provides the relevant built-in functionality to perform marking of an integer, via the $ReflectAboutInteger$ operation from the $Microsoft.Quantum.Arithmetic$. This means we no longer need to do it manually, but can instead just rely on the language feature instead.

Taking all of this into account, we can summarize our entire code to the following operation, which can be invoked from the Q# program the same way as we did it with the initial variant.

operation GroverVersion2(markIndex : Int, numberOfQubits : Int, iterations : Int) : Int {
    use qubits = Qubit[numberOfQubits];
    let register = LittleEndian(qubits);

    // superposition
    ApplyToEachA(H, qubits);
    DumpMachine();

    for x in 1..iterations {
        // mark the required number
        ReflectAboutInteger(markIndex, register);

        // amplitude amplification
        AmplifyAmplitude(qubits);
    }

    let number = MeasureInteger(register);
    ResetAll(qubits);

    return number;
}

The results produced upon the execution of this code should be identical to the previous version, with the obvious advantage that the code is now a lot more succinct.

Summary 🔗

In this post we generalized Grover’s algorithm and provided an implementation of it in Q# code. We tested the algorithm for some simple data sets and should be generally speaking pleased with the results, as they provide considerable speed up over the classical counterparts.

We can cap off our work with the quantum search with a very insightful quote from David Mermin’s book, which provides a bit of insight into possible use cases for this search method using today’s quantum hardware. Even though the problem define by Grover is often defined as “database search”, Mermin reminds us that:

Using as precious a resource as Qbits, however, merely to store classical information would be insanely extravagant, given our current or even our currently foreseeable ability to manufacture Qbits. Finding a unique solution – or one of a small number of solutions (…) – to a tough mathematical puzzle seems a more promising application.

This of course doesn’t render this technique useless, but highlights the fact, that with the current input/output limitations of quantum hardware, Grover’s search is best suited for the problem category where the data set can be artificially created in memory (e.g. a range of number).

We shall continue looking at quantum algorithms in the next part of this series!