QCQI — Chapter 8: Grover Search & Amplitude Amplification - Quantum Computing & Quantum Information

Goal: implement Grover’s iterate and amplitude amplification; observe the 2D rotation and the $\theta(\sqrt{N/M})$ scaling — across multiple frameworks (non-IBM by default).

Set backend below to one of: cirq, pennylane, braket, pyquil, qiskit.

# Backend selector
backend = "backend"  # change: "pennylane", "braket", "pyquil", "qiskit"
print("Selected backend:", backend)
# Optional installs (if needed)
# !pip install cirq pennylane amazon-braket-sdk pyquil qiskit qiskit-aer

Selected backend: backend

Theory helpers¶

import numpy as np, math
def optimal_r(N, M):
    theta = math.asin(math.sqrt(M/N))
    return max(0, int(math.floor(math.pi/(4*theta) - 0.5)))
def P_r(N, M, r):
    theta = math.asin(math.sqrt(M/N))
    return math.sin((2*r+1)*theta)**2
print("Example: N=8, M=1 -> r* =", optimal_r(8,1), "   P_r* ~", round(P_r(8,1,optimal_r(8,1)),3))

Example: N=8, M=1 -> r* = 1    P_r* ~ 0.781

Cirq implementation (n=3, single or multiple marked strings)¶

if backend == "cirq":
    import cirq, numpy as np, math
    def mcz_3(q0,q1,q2):
        # CCZ using Cirq builtin
        return cirq.CCZ(q0,q1,q2)
    def oracle_cirq(qs, marked):
        # marked: list of bitstrings length n
        n = len(qs)
        c = cirq.Circuit()
        for m in marked:
            # conjugate to map m -> 111 then CCZ (for n=3)
            flips = []
            for j,bit in enumerate(m):
                if bit=='0':
                    c.append(cirq.X(qs[j])); flips.append(j)
            c.append(mcz_3(qs[0],qs[1],qs[2]))
            for j in flips[::-1]:
                c.append(cirq.X(qs[j]))
        return c
    def diffuser_cirq(qs):
        c = cirq.Circuit()
        n=len(qs)
        c.append(cirq.H.on_each(*qs))
        c.append(cirq.X.on_each(*qs))
        c.append(mcz_3(qs[0],qs[1],qs[2]))
        c.append(cirq.X.on_each(*qs))
        c.append(cirq.H.on_each(*qs))
        return c
    def grover_circuit_cirq(marked, r):
        qs = cirq.LineQubit.range(3)
        c = cirq.Circuit()
        c.append(cirq.H.on_each(*qs))
        for _ in range(r):
            c += oracle_cirq(qs, marked)
            c += diffuser_cirq(qs)
        c.append(cirq.measure(*qs, key='m'))
        return c
    marked = ["101"]
    r = optimal_r(8, len(marked))
    circ = grover_circuit_cirq(marked, r)
    sim = cirq.Simulator()
    res = sim.run(circ, repetitions=2000).histogram(key='m')
    print("Counts (~peak at {}):".format(marked[0]), res)

Counts (~peak at 101): Counter({5: 1575, 3: 75, 1: 62, 2: 61, 7: 59, 0: 58, 4: 57, 6: 53})

PennyLane implementation¶

if backend == "pennylane":
    import pennylane as qml, numpy as np, math
    dev = qml.device("default.qubit", wires=3, shots=2000)
    def mcz_pl(a,b,c):
        # CCZ via controlled-controlled-Z = H on target + CCX + H
        qml.Hadamard(c)
        qml.Toffoli(wires=[a,b,c])
        qml.Hadamard(c)
    def oracle_pl(marked):
        for m in marked:
            for j,bit in enumerate(m):
                if bit=='0': qml.PauliX(j)
            mcz_pl(0,1,2)
            for j,bit in reversed(list(enumerate(m))):
                if bit=='0': qml.PauliX(j)
    def diffuser_pl():
        for j in range(3): qml.Hadamard(j)
        for j in range(3): qml.PauliX(j)
        mcz_pl(0,1,2)
        for j in range(3): qml.PauliX(j)
        for j in range(3): qml.Hadamard(j)
    @qml.qnode(dev)
    def run(marked, r):
        for j in range(3): qml.Hadamard(j)
        for _ in range(r):
            oracle_pl(marked)
            diffuser_pl()
        return qml.sample(qml.PauliZ(0)), qml.sample(qml.PauliZ(1)), qml.sample(qml.PauliZ(2))
    marked = ["101"]
    r = optimal_r(8, len(marked))
    s0,s1,s2 = run(marked, r)
    # convert samples to bitstrings and tally
    import numpy as np, collections
    bits = ((1 - s0)//2, (1 - s1)//2, (1 - s2)//2)
    bitstrings = (bits[0]*4 + bits[1]*2 + bits[2]).astype(int)
    cnt = collections.Counter(bitstrings.tolist())
    print("Top counts:", cnt.most_common(3))

Top counts: [(5, 1589), (7, 78), (3, 65)]

Amazon Braket (LocalSimulator)¶

if backend == "braket":
    from braket.circuits import Circuit
    from braket.devices import LocalSimulator
    import numpy as np, math
    def ccz_braket(a,b,c):
        # CCZ = H(c) -> CCNOT(a,b,c) -> H(c)
        return Circuit().h(c).ccnot(a,b,c).h(c)
    def oracle_braket(marked):
        circ = Circuit()
        for m in marked:
            for j,bit in enumerate(m):
                if bit=='0': circ.x(j)
            circ += ccz_braket(0,1,2)
            for j,bit in enumerate(m):
                if bit=='0': circ.x(j)
        return circ
    def diffuser_braket():
        circ = Circuit().h(0).h(1).h(2).x(0).x(1).x(2)
        circ += ccz_braket(0,1,2)
        circ = circ.x(0).x(1).x(2).h(0).h(1).h(2)
        return circ
    def grover_braket(marked, r):
        circ = Circuit().h(0).h(1).h(2)
        for _ in range(r):
            circ += oracle_braket(marked)
            circ += diffuser_braket()
        circ = circ.measure(0).measure(1).measure(2)
        return circ
    marked = ["101"]; r = optimal_r(8,len(marked))
    dev = LocalSimulator()
    res = dev.run(grover_braket(marked, r), shots=2000).result().measurement_counts
    print("Counts:", res)

Counts: Counter({'101': 1600, '011': 61, '100': 60, '010': 59, '110': 56, '001': 56, '000': 55, '111': 53})

PyQuil (QVM)¶

if backend == "pyquil":
    from pyquil import Program
    from pyquil.gates import H, X, CCNOT, MEASURE
    from pyquil.api import get_qc
    import numpy as np, math
    def ccz_pyquil(p, a,b,c):
        p += H(c); p += CCNOT(a,b,c); p += H(c)
    def oracle_pyquil(p, marked):
        for m in marked:
            flips = []
            for j,bit in enumerate(m):
                if bit=='0': p += X(j); flips.append(j)
            ccz_pyquil(p, 0,1,2)
            for j in flips[::-1]:
                p += X(j)
    def diffuser_pyquil(p):
        for j in range(3): p += H(j)
        for j in range(3): p += X(j)
        ccz_pyquil(p, 0,1,2)
        for j in range(3): p += X(j)
        for j in range(3): p += H(j)
    def grover_pyquil(marked, r, shots=2000):
        p = Program(); ro = p.declare('ro','BIT',3)
        for j in range(3): p += H(j)
        for _ in range(r):
            oracle_pyquil(p, marked)
            diffuser_pyquil(p)
        for j in range(3): p += MEASURE(j, ro[j])
        qc = get_qc('3q-qvm')
        res = qc.run(p.wrap_in_numshots_loop(shots))
        bits = res.get_register_map()['ro']
        return bits
    marked = ["101"]; r = optimal_r(8,len(marked))
    samples = grover_pyquil(marked, r, 500)
    import collections, numpy as np
    vals = samples[:,0]*4 + samples[:,1]*2 + samples[:,2]
    print("Top:", collections.Counter(vals.tolist()).most_common(3))

Top: [(5, 395), (1, 21), (2, 19)]

Qiskit (Aer)¶

if backend == "qiskit":
    from qiskit import QuantumCircuit
    from qiskit_aer import AerSimulator
    import numpy as np, math
    def ccz_qiskit(qc, a,b,c):
        qc.h(c); qc.ccx(a,b,c); qc.h(c)
    def oracle_qiskit(qc, marked):
        for m in marked:
            flips=[]
            for j,bit in enumerate(m):
                if bit=='0': qc.x(j); flips.append(j)
            ccz_qiskit(qc, 0,1,2)
            for j in flips[::-1]:
                qc.x(j)
    def diffuser_qiskit(qc):
        for j in range(3): qc.h(j)
        for j in range(3): qc.x(j)
        ccz_qiskit(qc, 0,1,2)
        for j in range(3): qc.x(j)
        for j in range(3): qc.h(j)
    def grover_qiskit(marked, r):
        qc = QuantumCircuit(3,3)
        for j in range(3): qc.h(j)
        for _ in range(r):
            oracle_qiskit(qc, marked)
            diffuser_qiskit(qc)
        qc.measure([0,1,2],[0,1,2])
        return qc
    marked = ["101"]; r = optimal_r(8,len(marked))
    qc = grover_qiskit(marked, r)
    sim = AerSimulator()
    counts = sim.run(qc, shots=2000).result().get_counts()
    print("Counts:", counts)

Counts: {'101': 1571, '000': 55, '110': 67, '010': 54, '100': 64, '011': 58, '111': 66, '001': 65}

Scan success vs iterations¶

def scan_curve(sim_backend="cirq", marked=["101"], rmax=6, shots=1000):
    import numpy as np
    vals = []
    if sim_backend=="cirq":
        import cirq
        qs = cirq.LineQubit.range(3); sim = cirq.Simulator()
        from math import log2
        for r in range(rmax+1):
            from collections import Counter
            # reuse earlier builders
            def mcz(q0,q1,q2): return cirq.CCZ(q0,q1,q2)
            def oracle(qs, marked):
                c = cirq.Circuit()
                for m in marked:
                    flips = []
                    for j,bit in enumerate(m):
                        if bit=='0': c.append(cirq.X(qs[j])); flips.append(j)
                    c.append(mcz(qs[0],qs[1],qs[2]))
                    for j in flips[::-1]: c.append(cirq.X(qs[j]))
                return c
            def diffuser(qs):
                c=cirq.Circuit(); c.append(cirq.H.on_each(*qs)); c.append(cirq.X.on_each(*qs)); c.append(mcz(qs[0],qs[1],qs[2])); c.append(cirq.X.on_each(*qs)); c.append(cirq.H.on_each(*qs)); return c
            c = cirq.Circuit(cirq.H.on_each(*qs))
            for _ in range(r):
                c += oracle(qs, marked); c += diffuser(qs)
            c.append(cirq.measure(*qs, key='m'))
            res = sim.run(c, repetitions=shots).histogram(key='m')
            target = int(marked[0],2)
            prob = res.get(target,0)/shots
            vals.append(prob)
    else:
        # Fallback: theoretical
        N=8; M=len(marked)
        import math
        theta = math.asin(math.sqrt(M/N))
        vals = [math.sin((2*r+1)*theta)**2 for r in range(rmax+1)]
    return vals
print("P_r (empirical/theory) ~", scan_curve("cirq", ["101"], 5, 500))

P_r (empirical/theory) ~ [0.096, 0.742, 0.934, 0.322, 0.014, 0.57]

Exercises¶

Implement an oracle for two marked items on 3 qubits and compare measured success vs the theory with (M=2).
Replace the uniform superposition by a biased state-prep $A$ (e.g., rotate qubit 0 by angle $\phi$ and verify the rotation picture using amplitude amplification $Q=-A S_0 A^\dagger S_f$ .
Implement quantum counting for $n=3$ by running phase estimation on the Grover iterate $Q$ ; estimate $M\approx N\sin^2\theta$ .