from functools import reduce
from typing import *
import time
import random
import math
import numpy as np
import matplotlib.pyplot as plt

from qiskit_aer import AerSimulator
from qiskit.quantum_info import Statevector, Operator
from qiskit.visualization import plot_histogram, plot_bloch_multivector
from qiskit import QuantumCircuit, transpile
sim = AerSimulator()

from util import zero, one, Z

QC: Grover’s Search

Under Construction

In this notebook, we’ll introduce Grover’s search. It enables us to search an unstructured dataset of length \(N\) in time \(O(\sqrt{N})\), which provides a quadratic speedup compared to a classical algorithm which takes \(O(N)\) time.

References

1. Introduction to Classical and Quantum Computing: Chapter 7.6

2. Qiskit notebook on Grover’s search

Problem: Unstructured Search

The problem of unstructured search involves searching a database, i.e., set of all bit strings \(\{0, 1\}^n\) over \(n\) bits for a particular bit string \(w\), called the marked bitstring. Intuitively, a classical algorithm would take \(O(N)\) time where \(N = 2^n\) since we must check every part of the database.

n = 3
N = 2**n
marked_bitstring = "010"
database = set([format(i, f"0{n}b") for i in range(2**n)])
database

{'000', '001', '010', '011', '100', '101', '110', '111'}

Reforumulation: Query

We can reformulate this as a query problem. In particular, define a function \(f_w: \{ 0, 1 \}^n \rightarrow \{0, 1\}\) such that \(f_w(x) = 1\) iff \(x = w\) and \(0\) otherwise. Then \(f_w\) can be put into one-to-one correspondence with databases \(\{0, 1\}^n\) with a marked string w.

def f_w(x) -> bool:
    return marked_string == x

Algorithm Walkthrough

We’ll now walkthrough Grover’s search algortihm.

marked_qubit = zero ^ one ^ zero
marked_qubit.draw("latex")

\[ |010\rangle\]

Step 1: Put first n qubits in superposition

Our first step is to put the \(n\) qubits in superposition. As usual, we will accomplish this by applying Hadamard gates to all qubits.

qc = QuantumCircuit(n, n)
qc.h(range(n))
qc.draw(output="mpl", style="iqp")

State after step 1

After the application of the H gates, we have the following state

\[ |s\rangle = H^{\otimes n} |0\rangle^{\otimes n} = |+\rangle^{\otimes n} \,. \]

s1 = (zero ^ zero ^ zero).evolve(Operator(qc))
s1.draw("latex")

\[\frac{\sqrt{2}}{4} |000\rangle+\frac{\sqrt{2}}{4} |001\rangle+\frac{\sqrt{2}}{4} |010\rangle+\frac{\sqrt{2}}{4} |011\rangle+\frac{\sqrt{2}}{4} |100\rangle+\frac{\sqrt{2}}{4} |101\rangle+\frac{\sqrt{2}}{4} |110\rangle+\frac{\sqrt{2}}{4} |111\rangle\]

State Representation

For Grover’s search, it is conveneient to consider a representation of the intermediate quantum states as a linear combination of \(|w\rangle\) and the states orthogonal to \(|w\rangle\). We can rewrite the state above as

\[\begin{align*} |s\rangle & = \frac{1}{\sqrt{N}} |w\rangle + \frac{1}{\sqrt{N}} \sum_{i \neq w} |i\rangle \tag{factor} \\ &= \frac{1}{\sqrt{N}} |w\rangle + \frac{\sqrt{N-1}}{\sqrt{N}} \underbrace{\frac{1}{\sqrt{N - 1}} \sum_{i \neq w} |i\rangle}_{|r\rangle} \tag{rearrange} \\ &= \underbrace{\frac{1}{\sqrt{N}}}_{\sin\theta} |w\rangle + \underbrace{\sqrt{\frac{N-1}{N}}}_{\cos \theta} |r\rangle \tag{substitute $|r\rangle$} \\ &= \sin \theta |w\rangle + \cos \theta |r\rangle \tag{$sin^2\theta + \cos^2 \theta = 1$} \end{align*}\]

where

1. \(|r\rangle = \frac{1}{\sqrt{N - 1}} \sum_{i \neq w} |i\rangle\),

2. \(\sin \theta = \frac{1}{\sqrt{N}}\), and

3. \(\cos \theta = \frac{\sqrt{N-1}}{\sqrt{N}}\).

Implicilty, \(|w\rangle\) uses the binary representation of \(w\), i.e., \(|w\rangle = |\overline{w}\rangle\). We do this for notational simplicity.

def mk_r(marked_qubit: Statevector) -> Statevector:
    N = len(marked_qubit)
    bitstrs = [format(i, f"0{n}b") for i in range(2**n)]
    comp_basis = [reduce(lambda x, y: x ^ y, [zero if x == "0" else one for x in b]) for b in bitstrs]
    without_marked = [b for b in comp_basis if not np.allclose(b, marked_qubit)]
    r = reduce(lambda x, y: x + y, without_marked)
    return 1/np.sqrt(N - 1) * r
mk_r(marked_qubit).draw("latex")

\[\frac{\sqrt{7}}{7} |000\rangle+\frac{\sqrt{7}}{7} |001\rangle+\frac{\sqrt{7}}{7} |011\rangle+\frac{\sqrt{7}}{7} |100\rangle+\frac{\sqrt{7}}{7} |101\rangle+\frac{\sqrt{7}}{7} |110\rangle+\frac{\sqrt{7}}{7} |111\rangle\]

((1/np.sqrt(N) * marked_qubit + np.sqrt((N - 1)/N)*mk_r(marked_qubit))).draw("latex")

\[\frac{\sqrt{2}}{4} |000\rangle+\frac{\sqrt{2}}{4} |001\rangle+\frac{\sqrt{2}}{4} |010\rangle+\frac{\sqrt{2}}{4} |011\rangle+\frac{\sqrt{2}}{4} |100\rangle+\frac{\sqrt{2}}{4} |101\rangle+\frac{\sqrt{2}}{4} |110\rangle+\frac{\sqrt{2}}{4} |111\rangle\]

Visualization

This state representation enables us to construct two helpful visualizations for seeing the quantum state as it evolve through Grover’s algorithm:

1. a angle visualization and

2. a amplitude and phase visualization.

Angle Visualization

We can visualize the angle between the current state \(|s\rangle\) and the marked state \(|w\rangle\). The smaller the angle, the more likely a measurement will collapse to the state \(|w\rangle\), our solution state.

def plot_angle(ax, marked_qubit: Statevector, s: Statevector) -> None:
    theta = np.pi/2 - np.arccos(np.dot(s, marked_qubit).real)
    print("cos(theta):", np.cos(theta), "sin(theta):", np.sin(theta))
    ax.set_xlim(-.1, 1)
    ax.set_ylim(-1.1, 1.1)
    ax.set_aspect("equal")
    ax.quiver(0, 0, np.cos(theta), np.sin(theta), angles='xy', scale_units='xy', scale=1, color='green', label=r'$|s\rangle$')
    ax.quiver(0, 0, 0, 1, angles='xy', scale_units='xy', scale=1, color='red', label=r'$|w\rangle$')
    ax.quiver(0, 0, 1, 0, angles='xy', scale_units='xy', scale=1, color='blue', label=r'$|r\rangle$')
    ax.legend()

fig, ax = plt.subplots(1, 1, figsize=(4, 4))
plot_angle(ax, marked_qubit, s1)

cos(theta): 0.9354143466934854 sin(theta): 0.3535533905932735

Amplitude and Phase Visualization

We introduce a visualization the tracks the “probability” amplitudes of the quantum state as Grover’s search is applied. The larger the magnitude of the probability amplitude, the more likely a measurement will collapse the current state to \(|w\rangle\), our solution state. The marked state that we are searching for is given in red. Every other state is given in blue.

def plot_amplitudes(ax, marked_state: str, s: Statevector) -> None:
    real_amps = [s.data[i].real for i in range(len(s))]
    N = len(real_amps)
    x_axis = np.arange(0, N, 1)
    colors = ["blue"] * N
    colors[int(marked_state, 2)] = "red"
    n = int(math.log2(N))
    axis_labels = [format(i, f"0{n}b") for i in range(2**n)]
    ax.grid()
    ax.bar(x_axis, real_amps, width=0.25, tick_label=axis_labels, color=colors)
    ax.set_ylim(-1, 1)

fig, ax = plt.subplots(1, 1, figsize=(8, 2))
plot_amplitudes(ax, marked_bitstring, s1)

Plotting together

We can put both plots side by side.

def plot_angles_and_amplitudes(marked_state: str, marked_qubit: Statevector, s: Statevector) -> None:
    fig, axs = plt.subplots(1, 2, figsize=(8, 4))
    plot_angle(axs[0], marked_qubit, s)
    plot_amplitudes(axs[1], marked_bitstring, s)

plot_angles_and_amplitudes(marked_bitstring, marked_qubit, s1)

cos(theta): 0.9354143466934854 sin(theta): 0.3535533905932735

Step 2: Iterative Amplification

We will iteratively apply the following two steps:

1. Grover’s oracle and

2. Grover’s diffuser.

Step 2.1: Apply Oracle

We define an oracle \(U_f\) that flips the phase of the marked state.

\[\begin{split} U_f |x\rangle = \begin{cases} -|w\rangle & \mbox{$x = w$} \\ -|x\rangle & \mbox{otherwise} \,. \end{cases} \end{split}\]

Towards this end, we will use the binary representation of \(|w\rangle\) as control bits to a circuit that will flip the phase

Aside: Multi-Controlled Z Gate (MCZ)

The \(MCZ\) gate is a gate that flips the target qubit’s phase if all the control qubits are set and does nothing otherwise. Unfortunately, Qiskit does not provide a MCZ gate so we will need to build one from scratch.

First, we can construct a Z Gate from X and H.

qc_z = QuantumCircuit(1)
qc_z.h(0)
qc_z.x(0)
qc_z.h(0)
qc_z.draw(output="mpl", style="iqp")

np.allclose(Operator(qc_z), Z)

True

Second, we do have access to a MCX gate, i.e., a multi-controlled X gate. This applies a X gate to a target qubit if all control qubits are set.

qc_mcx = QuantumCircuit(3)
qc_mcx.mcx([0, 1], 2)   # controls, target
qc_mcx.draw(output="mpl", style="iqp")

# MCX|011> = |111>
(zero ^ one ^ one).evolve(Operator(qc_mcx)).draw("latex")

\[ |111\rangle\]

Finally, we can put this together.

qc_mcz = QuantumCircuit(3)
qc_mcz.h(2)
qc_mcz.mcx([0, 1], 2)   # controls, target
qc_mcz.h(2)
qc_mcz.draw(output="mpl", style="iqp")

Operator(qc_mcz).draw("latex")

\[\begin{split}\begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & -1 \\ \end{bmatrix} \end{split}\]

# MCZ|111> = -|111>
(one ^ one ^ one).evolve(Operator(qc_mcz)).draw("latex")

\[- |111\rangle\]

Grover’s oracle from MCZ

The MCZ gate flipped the phase of the state \(|111\rangle\). In order to flip the phase of the marked state, we need to map it to the state \(|111\rangle\) and then “undo” our mapping which is possible since quantum operations are reversible. Towards this end, we observe that if a marked state contains a \(0\) we should apply an X gate to set it. For our target qubit, if it is \(0\) then the \(X\) gate will flip it to \(1\) and if it is \(1\), then it already has the correct bit set.

def grover_oracle(qc, marked_state):    
    n_qubits = qc.num_qubits
    # Apply X gates to qubits that should be |0⟩ in the marked state
    for i, bit in enumerate(marked_state):
        if bit == "0":
            qc.x(i)
    
    # ---- MCZ -----
    qc.h(n_qubits - 1)  # Change the target qubit to the |-| state
    qc.mcx(list(range(n_qubits - 1)), n_qubits - 1)  # Multi-controlled X with n-1 controls
    qc.h(n_qubits - 1)  # Undo the |-| state
    # ---- MCZ -----
    
    # Apply X gates again to revert the qubits to their original state
    for i, bit in enumerate(marked_state):
        if bit == '0':
            qc.x(i)

qc.barrier(label="initial | amp ")
grover_oracle(qc, marked_bitstring)
qc.draw(output="mpl", style="iqp")

s2_1 = (zero ^ zero ^ zero).evolve(Operator(qc))
s2_1.draw("latex")

\[\frac{\sqrt{2}}{4} |000\rangle+\frac{\sqrt{2}}{4} |001\rangle- \frac{\sqrt{2}}{4} |010\rangle+\frac{\sqrt{2}}{4} |011\rangle+\frac{\sqrt{2}}{4} |100\rangle+\frac{\sqrt{2}}{4} |101\rangle+\frac{\sqrt{2}}{4} |110\rangle+\frac{\sqrt{2}}{4} |111\rangle\]

plot_angles_and_amplitudes(marked_bitstring, marked_qubit, s2_1)

cos(theta): 0.9354143466934853 sin(theta): -0.35355339059327373

Step 2.2: Diffuser

We now would like to boost the probability of the marked state somehow. This should be possible since we have flipped its phase so it is somehow distinguished from the other states. We accomplish this with the following operator

\[ R_w = 2|w\rangle\langle w| - I \,. \]

This leaves \(|w\rangle\) alone

\[\begin{align*} R_w |w\rangle & = (2|w\rangle\langle w| - I)|w\rangle \\ & = 2|w\rangle\langle w|w\rangle - |w\rangle \\ & = 2|w\rangle - |w\rangle \\ & = |w\rangle \,. \end{align*}\]

Additionally, it flips the phase of any state vector orthogonal to \(|w\rangle\) as in

\[\begin{align*} R_w |w^\perp\rangle & = (2|w\rangle\langle w| - I)||w^\perp\rangle \\ & = 2|w\rangle\langle w|w^\perp\rangle - |w^\perp\rangle \\ & = -|w^\perp\rangle \,. \end{align*}\]

Circuit for the diffuser

Observe that

\[\begin{align*} R_w = H^{\otimes n} \underbrace{(2|0\rangle\langle 0| - I)}_{R_0} H^{\otimes n} \,. \end{align*}\]

Thus the problem reduces to constructing \(R_0\).

def grover_diffuser(qc):    
    n_qubits = qc.num_qubits
    
    # H^n
    qc.h(range(n_qubits))

    # R_0
    qc.x(range(n_qubits))
    qc.h(n_qubits - 1)
    qc.mcx(list(range(n_qubits - 1)), n_qubits - 1)
    qc.h(n_qubits - 1)
    qc.x(range(n_qubits))
    
    # H^n
    qc.h(range(n_qubits))
    
    return qc

grover_diffuser(qc)
qc.draw(output="mpl", style="iqp")

s2_2 = (zero ^ zero ^ zero).evolve(Operator(qc))
s2_2.draw("latex")

\[-0.1767766953 |000\rangle-0.1767766953 |001\rangle-0.8838834765 |010\rangle-0.1767766953 |011\rangle-0.1767766953 |100\rangle-0.1767766953 |101\rangle-0.1767766953 |110\rangle-0.1767766953 |111\rangle\]

plot_angles_and_amplitudes(marked_bitstring, marked_qubit, s2_2)

cos(theta): 0.4677071733467445 sin(theta): -0.8838834764831834

How many iterations?

We see geometrically that each iteration of applying Grover’s oracle and diffuser increases the angle between \(|s\rangle\) and \(|r\rangle\) by \(\theta\). Consequently, after iteration, we should get closer and closer to the solution state. After some algebra, we see that the optimal number of iterations \(T\) is obtained by solving

\[\begin{align*} \theta + 2T\theta = \frac{\pi}{2} \,. \end{align*}\]

Thus the optimal number of iteration is approximately $\( T \approx \frac{\pi}{4}\sqrt{N} \,. \)$

qc.barrier(label=f"End Iter {1}")
grover_oracle(qc, marked_bitstring)
qc.draw(output="mpl", style="iqp")

s2_iter2_1 = (zero ^ zero ^ zero).evolve(Operator(qc))
plot_angles_and_amplitudes(marked_bitstring, marked_qubit, s2_iter2_1)

cos(theta): 0.4677071733467449 sin(theta): 0.8838834764831832

grover_diffuser(qc)
qc.barrier(label=f"End Iter {2}")
qc.draw(output="mpl", style="iqp")

s2_iter2_2 = (zero ^ zero ^ zero).evolve(Operator(qc))
s2_iter2_2.draw("latex")

\[-0.0883883476 |000\rangle-0.0883883476 |001\rangle+0.9722718241 |010\rangle-0.0883883476 |011\rangle-0.0883883476 |100\rangle-0.0883883476 |101\rangle-0.0883883476 |110\rangle-0.0883883476 |111\rangle\]

plot_angles_and_amplitudes(marked_bitstring, marked_qubit, s2_iter2_2)

cos(theta): 0.23385358667338038 sin(theta): 0.9722718241315007

Step 3: Measure

Observe that the amplitude of the marked state \(|010\rangle\) is \(0.97\). This means the probability of measure it is about 95%. Thus the last step of Grover’s is to measure and extract the state.

def grover_measure(qc):
    for n in range(qc.num_qubits):
        qc.measure(n, n)

grover_measure(qc)
qc.draw(output="mpl", style="iqp")

t_qc = transpile(qc, sim)
results = sim.run(t_qc, shots=1024).result()
counts = results.get_counts()
plot_histogram(counts)

Putting it all together

We can put this all together now.

def grovers(n_qubit: int, marked_bitstring: str) -> QuantumCircuit:
    qc = QuantumCircuit(n_qubit, n_qubit)
    
    # Step 1: Superposition
    qc.h(range(n_qubit))
    
    # Step 2: Iterative Amplification
    num_iter = int(np.floor(np.pi / 4 * np.sqrt(2 ** n_qubit)))
    qc.barrier(label="Step 1|Step 2: Diffuser")
    for i in range(num_iter):
        # Step 2.1: Oracle
        grover_oracle(qc, marked_bitstring)
        qc.barrier(label=f"Oracle | Diffuser")
        
        # Step 2.2: Diffuser
        grover_diffuser(qc)
        qc.barrier(label=f"End Iter {i+1}")

    # Step 3: Read out results
    grover_measure(qc)
        
    return qc

qc = grovers(3, "111")
qc.draw(output="mpl", style="iqp")

t_qc = transpile(qc, sim)
results = sim.run(t_qc, shots=1024).result()
counts = results.get_counts()
plot_histogram(counts)

Summary

1. We introduced Grover’s search, a quantum algorithm for unstructured search.

2. Grover’s search provides a quadratic speedup compared to classical search.