Grover Search Algorithm

This notebook is an adapted version from https://github.com/QISKit/qiskit-tutorial. We show to perform the Grover Search algorithm both on a local simulator and on the Quantum Inspire backend.

For more information about how to use the IBM Q Experience (QX), consult the tutorials, or check out the community.

Contributors Pieter Eendebak, Giacomo Nannicini and Rudy Raymond (based on this article)

Introduction

Grover search is one of the most popular algorithms used for searching a solution among many possible candidates using Quantum Computers. If there are $N$ possible solutions among which there is exactly one solution (that can be verified by some function evaluation), then Grover search can be used to find the solution with $O(\sqrt{N})$ function evaluations. This is in contrast to classical computers that require $\Omega(N)$ function evaluations: the Grover search is a quantum algorithm that provably can be used search the correct solutions quadratically faster than its classical counterparts.

Here, we are going to illustrate the use of Grover search to find a particular value in a binary number. The key elements of Grovers algorithm are:

1. Initialization to a uniform superposition
2. The oracle function
3. Reflections (amplitude amplification)

import numpy as np
import os

from qiskit.tools.visualization import plot_histogram
from qiskit import execute, QuantumCircuit, QuantumRegister, ClassicalRegister
from qiskit import BasicAer
from IPython.display import display, Math

from quantuminspire.credentials import get_authentication
from quantuminspire.qiskit import QI

QI_URL = os.getenv('API_URL', 'https://api.quantum-inspire.com/')

The oracle function

We implement an oracle function (black box) that acts as -1 on a single basis state, and +1 on all other status.

def format_vector(state_vector, decimal_precision=7):
    """ Format the state vector into a LaTeX formatted string.

    Args:
        state_vector (list or array): The state vector with complex
                                      values e.g. [-1, 2j+1].

    Returns:
        str: The LaTeX format.
    """
    result = []
    epsilon = 1/pow(10, decimal_precision)
    bit_length = (len(state_vector) - 1).bit_length()
    for index, complex_value in enumerate(state_vector):
        has_imag_part = np.round(complex_value.imag, decimal_precision) != 0.0
        value = complex_value if has_imag_part else complex_value.real
        value_round = np.round(value, decimal_precision)
        if np.abs(value_round) < epsilon:
            continue

        binary_state = '{0:0{1}b}'.format(index, bit_length)
        result.append(r'{0:+2g}\left\lvert {1}\right\rangle '.format(value_round, binary_state))
    return ''.join(result)

def run_circuit(q_circuit, q_register, number_of_qubits=None, backend_name='statevector_simulator'):
    """ Run a circuit on all base state vectors and show the output.

    Args:
        q_circuit (QuantumCircuit):
        q_register (QuantumRegister)
        number_of_qubits (int or None): The number of qubits.
        backend (str): ...
    """
    if not isinstance(number_of_qubits, int):
        number_of_qubits = q_register.size

    if q_register.size != number_of_qubits:
        warnings.warn('incorrect register size?')

    latex_text = r'\mathrm{running\ circuit\ on\ set\ of\ basis\ states:}'
    display(Math(latex_text))

    base_states = 2 ** number_of_qubits
    backend = BasicAer.get_backend(backend_name)
    for base_state in range(base_states):
        pre_circuit = QuantumCircuit(q_register)
        state = base_state
        for kk in range(number_of_qubits):
            if state % 2 == 1:
                pre_circuit.x(q[kk])
            state = state // 2

        input_state = r'\left\lvert{0:0{1}b}\right\rangle'.format(base_state, number_of_qubits)
        circuit_total = pre_circuit.compose(q_circuit)
        job = execute(circuit_total, backend=backend)
        output_state = job.result().get_statevector(circuit_total)

        latex_text = input_state + r'\mathrm{transforms\ to}: ' + format_vector(output_state)
        display(Math(latex_text))

n=3
N=2**n
q = QuantumRegister(n)
qc = QuantumCircuit(q)

if n==1:
    def black_box(qc, q):
        qc.z(q)
elif n==2:
    def black_box(qc, q):
        for i in range(n):
            qc.s(q[i])
        qc.h(q[1])
        qc.cx(q[0], q[1])
        qc.h(q[1])
        for i in range(n):
            qc.s(q[i])
else:
    def black_box(qc, q):
        qc.h(q[2])
        qc.ccx(q[0], q[1], q[2])
        qc.h(q[2])
black_box(qc,q)
cplot=qc.draw(output='mpl')
display(cplot)

print('black box circuit:')
run_circuit(qc, q)

black box circuit:

$\displaystyle \mathrm{running\ circuit\ on\ set\ of\ basis\ states:}$

$\displaystyle \left\lvert000\right\rangle\mathrm{transforms\ to}: +1\left\lvert 000\right\rangle $

$\displaystyle \left\lvert001\right\rangle\mathrm{transforms\ to}: +1\left\lvert 001\right\rangle $

$\displaystyle \left\lvert010\right\rangle\mathrm{transforms\ to}: +1\left\lvert 010\right\rangle $

$\displaystyle \left\lvert011\right\rangle\mathrm{transforms\ to}: +1\left\lvert 011\right\rangle $

$\displaystyle \left\lvert100\right\rangle\mathrm{transforms\ to}: +1\left\lvert 100\right\rangle $

$\displaystyle \left\lvert101\right\rangle\mathrm{transforms\ to}: +1\left\lvert 101\right\rangle $

$\displaystyle \left\lvert110\right\rangle\mathrm{transforms\ to}: +1\left\lvert 110\right\rangle $

$\displaystyle \left\lvert111\right\rangle\mathrm{transforms\ to}: -1\left\lvert 111\right\rangle $

Inversion about the average

Another important procedure in Grover search is to have an operation that perfom the inversion-about-the-average step, namely, it performs the following transformation:

$$ \sum_{j=0}^{2^{n}-1} \alpha_j |j\rangle \rightarrow \sum_{j=0}^{2^{n}-1}\left(2 \left( \sum_{k=0}^{k=2^{n}-1} \frac{\alpha_k}{2^n} \right) - \alpha_j \right) |j\rangle $$

The above transformation can be used to amplify the probability amplitude $\alpha_s$ when s is the solution and $\alpha_s$ is negative (and small), while $\alpha_j$ for $j \neq s$ is positive. Roughly speaking, the value of $\alpha_s$ increases by twice the average of the amplitudes, while others are reduced. The inversion-about-the-average can be realized with the sequence of unitary matrices as below:

$$ H^{\otimes n} \left(2|0\rangle \langle 0 | - I \right) H^{\otimes n} $$

The first and last $H$ are just Hadamard gates applied to each qubit. The operation in the middle requires us to design a sub-circuit that flips the probability amplitude of the component of the quantum state corresponding to the all-zero binary string. The sub-circuit can be realized by the following function, which is a multi-qubit controlled-Z which flips the probability amplitude of the component of the quantum state corresponding to the all-one binary string. Applying X gates to all qubits before and after the function realizes the sub-circuit.

def n_controlled_Z(circuit, controls, target):
    """Implement a Z gate with multiple controls"""
    if (len(controls) > 2):
        raise ValueError('The controlled Z with more than 2 ' +
                         'controls is not implemented')
    elif (len(controls) == 1):
        circuit.h(target)
        circuit.cx(controls[0], target)
        circuit.h(target)
    elif (len(controls) == 2):
        circuit.h(target)
        circuit.ccx(controls[0], controls[1], target)
        circuit.h(target)

Finally, the inversion-about-the-average circuit can be realized by the following function:

def inversion_about_average(circuit, f_in, n):
    """Apply inversion about the average step of Grover's algorithm."""
    # Hadamards everywhere
    if n==1:
        circuit.x(f_in[0])
        return
    for j in range(n):
        circuit.h(f_in[j])
    # D matrix: flips the sign of the state |000> only
    for j in range(n):
        circuit.x(f_in[j])
    n_controlled_Z(circuit, [f_in[j] for j in range(n-1)], f_in[n-1])
    for j in range(n):
        circuit.x(f_in[j])
    # Hadamards everywhere again
    for j in range(n):
        circuit.h(f_in[j])

We show the circuit that performs inversion about the average on $n$ qubits. We also show the effect of the circuit on the basis states.

qInvAvg = QuantumCircuit(q)
inversion_about_average(qInvAvg, q, n)
qInvAvg.draw(output='mpl')

print('inversion average circuit:')
qc = QuantumCircuit(q)
inversion_about_average(qc, q, n)
run_circuit(qc, q, n)

inversion average circuit:

$\displaystyle \mathrm{running\ circuit\ on\ set\ of\ basis\ states:}$

$\displaystyle \left\lvert000\right\rangle\mathrm{transforms\ to}: +0.75\left\lvert 000\right\rangle -0.25\left\lvert 001\right\rangle -0.25\left\lvert 010\right\rangle -0.25\left\lvert 011\right\rangle -0.25\left\lvert 100\right\rangle -0.25\left\lvert 101\right\rangle -0.25\left\lvert 110\right\rangle -0.25\left\lvert 111\right\rangle $

$\displaystyle \left\lvert001\right\rangle\mathrm{transforms\ to}: -0.25\left\lvert 000\right\rangle +0.75\left\lvert 001\right\rangle -0.25\left\lvert 010\right\rangle -0.25\left\lvert 011\right\rangle -0.25\left\lvert 100\right\rangle -0.25\left\lvert 101\right\rangle -0.25\left\lvert 110\right\rangle -0.25\left\lvert 111\right\rangle $

$\displaystyle \left\lvert010\right\rangle\mathrm{transforms\ to}: -0.25\left\lvert 000\right\rangle -0.25\left\lvert 001\right\rangle +0.75\left\lvert 010\right\rangle -0.25\left\lvert 011\right\rangle -0.25\left\lvert 100\right\rangle -0.25\left\lvert 101\right\rangle -0.25\left\lvert 110\right\rangle -0.25\left\lvert 111\right\rangle $

$\displaystyle \left\lvert011\right\rangle\mathrm{transforms\ to}: -0.25\left\lvert 000\right\rangle -0.25\left\lvert 001\right\rangle -0.25\left\lvert 010\right\rangle +0.75\left\lvert 011\right\rangle -0.25\left\lvert 100\right\rangle -0.25\left\lvert 101\right\rangle -0.25\left\lvert 110\right\rangle -0.25\left\lvert 111\right\rangle $

$\displaystyle \left\lvert100\right\rangle\mathrm{transforms\ to}: -0.25\left\lvert 000\right\rangle -0.25\left\lvert 001\right\rangle -0.25\left\lvert 010\right\rangle -0.25\left\lvert 011\right\rangle +0.75\left\lvert 100\right\rangle -0.25\left\lvert 101\right\rangle -0.25\left\lvert 110\right\rangle -0.25\left\lvert 111\right\rangle $

$\displaystyle \left\lvert101\right\rangle\mathrm{transforms\ to}: -0.25\left\lvert 000\right\rangle -0.25\left\lvert 001\right\rangle -0.25\left\lvert 010\right\rangle -0.25\left\lvert 011\right\rangle -0.25\left\lvert 100\right\rangle +0.75\left\lvert 101\right\rangle -0.25\left\lvert 110\right\rangle -0.25\left\lvert 111\right\rangle $

$\displaystyle \left\lvert110\right\rangle\mathrm{transforms\ to}: -0.25\left\lvert 000\right\rangle -0.25\left\lvert 001\right\rangle -0.25\left\lvert 010\right\rangle -0.25\left\lvert 011\right\rangle -0.25\left\lvert 100\right\rangle -0.25\left\lvert 101\right\rangle +0.75\left\lvert 110\right\rangle -0.25\left\lvert 111\right\rangle $

$\displaystyle \left\lvert111\right\rangle\mathrm{transforms\ to}: -0.25\left\lvert 000\right\rangle -0.25\left\lvert 001\right\rangle -0.25\left\lvert 010\right\rangle -0.25\left\lvert 011\right\rangle -0.25\left\lvert 100\right\rangle -0.25\left\lvert 101\right\rangle -0.25\left\lvert 110\right\rangle +0.75\left\lvert 111\right\rangle $

Grover Search: putting all together

The complete steps of Grover search is as follow.

1. Create the superposition of all possible solutions as the initial state (with working qubits initialized to zero) $$ \sum_{j=0}^{2^{n}-1} \frac{1}{2^n} |j\rangle |0\rangle$$

2. Repeat for $T$ times:

   • Apply the blackbox function
   • Apply the inversion-about-the-average function

3. Measure to obtain the solution

Before we go to the code to perform the Grover search we make some remarks on the number of repetitions $T$ that we have to perform (for details see Grover algorithm, Wikipedia).

Each Grover step rotates the 'winner solution' by a fixed angle. This means that after a certain number of steps we arrive at the optimal approximation (e.g. the amplitude of the winner solution is maximal). If we then apply more iterations, the quality of our result will go down. For a database of size $N=2^n$ the optimal number of iterations is $$r=\pi \sqrt{N}/4$$

theta = 2*np.arcsin(1/np.sqrt(N))
r=np.pi*np.sqrt(N)/4
display(Math(r'\textrm{Rotation of the winner: } \theta = %.2f \mathrm{\ [deg]}' % (np.rad2deg(theta))) )
print('Optimal number of Grover iterations for n=%d: %.1f' % (n,r) )
T=int(r)

$\displaystyle \textrm{Rotation of the winner: } \theta = 41.41 \mathrm{\ [deg]}$

Optimal number of Grover iterations for n=3: 2.2

The probablity of the winner state after $T$ iterations is $\sin( (T+1/2)\theta)^2$

for i in range(int(r+2)):
    p=np.sin((i+1/2)*theta)**2
    print('%d iterations: p %.2f' % (i, p))

0 iterations: p 0.12
1 iterations: p 0.78
2 iterations: p 0.95
3 iterations: p 0.33

Finally we define the complete circuit for Grovers algorithm, excute it and show the results.

"""Grover search implemented in QISKit.

This module contains the code necessary to run Grover search on 3
qubits, both with a simulator and with a real quantum computing
device. This code is the companion for the paper
"An introduction to quantum computing, without the physics",
Giacomo Nannicini, https://arxiv.org/abs/1708.03684.

"""
def input_state(circuit, f_in, n):
    """(n+1)-qubit input state for Grover search."""
    for j in range(n):
        circuit.h(f_in[j])

q = QuantumRegister(n)
ans = ClassicalRegister(n)
qc = QuantumCircuit(q, ans)

input_state(qc, q, n)

backend=BasicAer.get_backend('statevector_simulator')
job = execute(qc, backend=backend, shots=10)
result = job.result()
state_vector = result.get_statevector(qc)
m=display( Math('\mathrm{state\ after\ initialization:\ }' +format_vector(state_vector)))

# apply T rounds of oracle and inversion about the average
print('number of iterations T=%d'% T)
for t in range(T):
    for i in range(n):
        qc.barrier(q[i]) # for better visualization
    qc.i(q[0])
    # Apply T full iterations
    black_box(qc, q)
    for i in range(n):
        qc.barrier(q[i])
    qc.i(q[0])
    inversion_about_average(qc, q, n)

# Measure the output register in the computational basis
for j in range(n):
    qc.measure(q[j], ans[j])

# Execute circuit
backend=BasicAer.get_backend('qasm_simulator')
job  = execute(qc, backend=backend, shots=10)
result = job.result()

# Get counts and plot histogram
counts = result.get_counts()
plot_histogram(counts)

$\displaystyle \mathrm{state\ after\ initialization:\ }+0.353553\left\lvert 000\right\rangle +0.353553\left\lvert 001\right\rangle +0.353553\left\lvert 010\right\rangle +0.353553\left\lvert 011\right\rangle +0.353553\left\lvert 100\right\rangle +0.353553\left\lvert 101\right\rangle +0.353553\left\lvert 110\right\rangle +0.353553\left\lvert 111\right\rangle $

number of iterations T=2

As expected, the state that is indicated by the oracle function has the highest probability of begin measured.

We show the full circuit that was generated by the code.

qc.draw(output='mpl')

Run the circuit on the Quantum Inspire simulator

First we make a connection to the Quantum Inspire website.

authentication = get_authentication()
QI.set_authentication(authentication, QI_URL)

We can list backends and perform other functions with the QuantumInspireProvider.

QI.backends()

[<QuantumInspireBackend('QX-34-L') from QI()>,
 <QuantumInspireBackend('Spin-2') from QI()>,
 <QuantumInspireBackend('QX single-node simulator') from QI()>,
 <QuantumInspireBackend('Starmon-5') from QI()>]

We create a QisKit backend for the Quantum Inspire interface and execute the circuit generated above.

qi_backend = QI.get_backend('QX single-node simulator')
j=execute(qc, backend=backend, shots=512)

We can wait for the results and then print them

result = j.result()
print('Generated histogram:')
print(result.get_counts())

Generated histogram:
{'111': 479, '010': 5, '100': 5, '110': 7, '001': 4, '101': 6, '000': 3, '011': 3}

Visualization can be done with the normal Python plotting routines, or with the QisKit SDK.

plot_histogram(result.get_counts(qc))

A screenshot from the execution result on the Quantum Inspire website.

References

[1] "A fast quantum mechanical algorithm for database search", L. K. Grover, Proceedings of the 28th Annual ACM Symposium on the Theory of Computing (STOC 1996)

[2] "Tight bounds on quantum searching", Boyer et al., Fortsch.Phys.46:493-506,1998

[3] "Quantum Inspire"