Grover’s Search

In this demonstration, we simulate Grover’s Search to find an element in an unstructured search space. We limit ourselves to use only the Pauli $X$, Hadamard and multi-controlled Pauli $Z$ gates.

Introduction #

Let $f(x)$ be a function which returns $0$ for every input $x$, except at value $x=\omega$, for which $f(\omega)=1$. We seek to find $\omega$ among the $N$ possible values of $x$, by querying $f(x)$. In a classical algorithm, we would have to try every value of $x$ in some arbitrary order until we discover $f(x) = 1$, hence the classical complexity is $O(N)$.

But assume we have access to an error-free quantum computer, and a unitary oracle $U_\omega$ which can flip the phase of our sought unknown state.
$$ \begin{cases} \; U_\omega |\omega \rangle = – |\omega \rangle & \\ \; U_\omega |x\rangle = \hphantom{-}|x\rangle & x \ne \omega \end{cases} $$ Armed with this utility, we can discover $\omega$ with a complexity of only $O(\sqrt{N})$ oracle queries, and $O(\log(N)\sqrt{N})$ quantum gates, using Grover’s algorithm.

Let $N = 2^{\text{numQubits}}$, implying any one of the possible computational basis states of our quantum computer could be the solution state $|\omega \rangle$. Grover’s algorithm proceeds as follows.

Implementation #

Let’s simulate it!

We first prepare the QuEST environment and our simulation parameters.

#include "QuEST.h"

int main() {
    
    // prepare the hardware-agnostic QuEST environment
    QuESTEnv env = createQuESTEnv()

    // choose the system size
    int numQubits = 15;
    int numElems = (int) pow(2, numQubits);

    // randomly choose the element for which to search
    srand(time(NULL));
    int omega = rand() % numElems;

    // prepare our register
    Qureg qureg = createQureg(numQubits, env);

    // search for omega using Grover's algorithm
    applyGrovers(qureg, omega);

    // tidy up
    destroyQureg(qureg, env);
    destroyQuESTEnv(env);
    return 0;
}

Algorithm #

With that out of the way, Grover’s algorithm is very simple.

void applyGrovers(Qureg qureg, int omega) {

    int numQubits = qureg.numQubitsRepresented;
    int numElems = (int) pow(2, numQubits);
    int numReps = ceil(M_PI/4 * sqrt(numElems));

    // prepare |+><+|
    initPlusState(qureg);

    // repeatedly apply Grover's iteration
    for (int r=0; r<numReps; r++) {
        applyOracle(qureg, numQubits, omega);
        applyDiffuser(qureg, numQubits);
    }

    // collapse into the output classical state
    measureResult(qureg, omega);
}

Oracle #

Our next task is to implement the oracle operator, which effects $$ \begin{cases} \; U_\omega |\omega \rangle = – |\omega \rangle & \\ \; U_\omega |x\rangle = \hphantom{-}|x\rangle & x \ne \omega \end{cases} $$ We could do this in a variety of ways, for example, by simply constructing a diagonal matrix which is almost identity, but flips the sign of $|\omega\rangle$
$$ U_\omega = \begin{pmatrix} 1 \\ & 1 \\ & & \ddots \\ & & & -1 \\ & & & & \ddots \\ & & & & & 1 \end{pmatrix} $$
Applying this via applyDiagonalOp would work, but would wastefully double our memory requirements.

// wasteful

int main() {

    // op is just as large as qureg!
    DiagonalOp op = createDiagonalOp(numQubits, env);
    for (int x=0; x<numElems; x++)
        op.real[x] = (x == omega)? -1 : 1;

    ...
}

void applyOracle(...) {
    applyDiagonalOp(qureg, op);
}

Another method would be to simply manually flip the sign of the single amplitude modified by our oracle.

void applyOracle(...) {

    Complex amp = getAmp(qureg, omega);
    amp.real *= -1;
    amp.imag *= -1;
    setAmps(qureg, omega, &(amp.real), &(amp.imag), 1);
}

This works just fine, though a pedant may wish to effect the oracle with elementary unitaries, even though this requires baking-in the solution. To do so, we first appreciate that the multi-controlled phase flip (a.k.a. controlled Pauli $Z$), with matrix form
$$ cZ = \begin{pmatrix} 1 \\ & 1 \\ & & \ddots \\ & & & 1 \\ & & & & -1 \end{pmatrix} $$ is similar to our $U_{\omega}$ oracle, but it flips the sign of the $|1\dots 1\rangle$ state, instead of $|\omega\rangle$. So, we simply first transform between states $|1\dots 1\rangle \; \leftrightarrow \; |\omega\rangle$ by flipping the zero-bits of $\omega$’s binary representation. We can then flip the phase of $|1\dots 1\rangle$ before transforming back.

/* effect |omega> -> -|omega> via a 
 * multi-controlled phase flip gate 
 */
void applyOracle(Qureg qureg, int numQubits, int omega) {
    
    // apply X to transform |111> into |omega>
    for (int q=0; q<numQubits; q++)
        if (((omega >> q) & 1) == 0) // q-th bit of omega
            pauliX(qureg, q);
        
    // effect |111> -> -|111>    
    int ctrls[numQubits];
    for (int q=0; q<numQubits; q++)
        ctrls[q] = q;
    multiControlledPhaseFlip(qureg, ctrls, numQubits);
    
    // apply X to transform |omega> into |111>
    for (int q=0; q<numQubits; q++)
        if (((omega >> q) & 1) == 0)  // q-th bit of omega
            pauliX(qureg, q);
}

Diffuser #

Now we implement the diffuser operator $$ 2\;|+\rangle\langle + | \; – \; I $$ There are again a multitude of ways we can effect this operator; we’ll make use of $X$, $H$ and $cZ$. We first observe that by transforming every qubit into the Hadamard basis, the diffuser becomes $$ 2\;|0\dots 0\rangle\langle 0 \dots 0| \; – \; I$$ This is again similar to the multi-controlled phase flip, which we can express as $$ c\dots cZ \;=\; \text{diag}(\;1, \; \dots, \; 1, \; -1 \; ) \;=\; I \; – \; 2 \; |1\dots 1\rangle\langle 1 \dots 1| $$ except that states $|0\dots 0\rangle \; \leftrightarrow \; |1\dots 1\rangle$ are swapped, and the global phase is flipped. Transforming between those states simply requires simply flipping all qubits with Pauli $X$.
$$ I \; – \; 2\;|0\dots 0\rangle\langle 0 \dots 0| \;\; = \;\; X_0 X_1 \dots \; ( c\dots cZ) \; \dots X_1 X_0 $$ This differs from the Hadamard basis diffuser only by an unphysical global phase of $e^{i \pi } = -1 $.

void applyDiffuser(Qureg qureg, int numQubits) {
    
    // apply H to transform |+> into |0>
    for (int q=0; q<numQubits; q++)
        hadamard(qureg, q);

    // apply X to transform |1..1> into |0..0>
    for (int q=0; q<numQubits; q++)
        pauliX(qureg, q);
    
    // effect |1..1> -> -|1..1>
    int ctrls[numQubits];
    for (int q=0; q<numQubits; q++)
        ctrls[q] = q;
    multiControlledPhaseFlip(qureg, ctrls, numQubits);
    
    // apply X to transform |0..0> into |1..1>
    for (int q=0; q<numQubits; q++)
        pauliX(qureg, q);
    
    // apply H to transform |0> into |+>
    for (int q=0; q<numQubits; q++)
        hadamard(qureg, q);
}

Naturally the adjacent $X$ and $H$ gates could be composed into a single gate $$ X \cdot H = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ -1 & 1 \end{pmatrix} $$ and effected via unitary or compactUnitary

qreal re = 1/sqrt(2);
Complex a = {.real =  re, .imag = 0};
Complex b = {.real = -re, .imag = 0};

for (int q=0; q<numQubits; q++)
    compactUnitary(qureg, q, a, b);

Measurement #

While repeatedly applying the oracle and diffuser operators, the quantum register should be converging into a single quantum state. At the end of the algorithm, an experimentalist would perform a final measurement on every qubit to infer the bits of $\omega$

void measureResult(...) {

    int omega = 0;

    for (int q=0; q<numQubits; q++) {

        int bit = measure(qureg, q);
        omega += bit * (int) pow(2, q);
    }
}

However, since $|\omega\rangle$ is a known computational basis state in our simulation, we can instead directly measure the probability of these measurements all collapsing into their correct states.

void measureResult(Qureg qureg, int omega) {

    qreal prob = getProbAmp(qureg, omega);

    printf("success probability: %g \n", prob);
}

If we did everything right, we should find that for $N \gg 1$, this probability converges to $1$. Grover’s algorithm works!

Run #

View the full code demonstration here.

Run the demonstration, from within the root QuEST directory, by running in terminal

cp examples/grovers_search.c .
make SOURCES=grovers_search
./demo