Chapter 33: Capstone: Your Quantum Algorithm Portfolio — Implemented, Executed on Real Hardware, and Analyzed

Learning Objectives

By the end of this capstone chapter, you will have:

• Implemented eight foundational quantum algorithms in Qiskit with complete, runnable code.
• Executed each algorithm on both a simulator and real IBM Quantum hardware.
• Analyzed the effects of noise, decoherence, and gate errors on real quantum results.
• Compared ideal (simulated) outcomes with noisy (hardware) outcomes.
• Built a portfolio of quantum programs demonstrating practical quantum computing skills.
• Understood the gap between theoretical quantum algorithms and their NISQ-era implementations.
• Computed quantitative noise budgets and fidelity estimates for each algorithm.
• Applied error mitigation techniques and measured their effectiveness.
• Developed intuition for which algorithms are NISQ-feasible and which require fault tolerance.

33.1 The Capstone Philosophy

This chapter is not a passive reading experience. It is a guided workshop. Each section presents a quantum algorithm, derives its mathematical foundation, provides a complete Qiskit implementation, and instructs you to run it on real IBM Quantum hardware. By the end, you will have executed the major quantum algorithms on actual quantum processors — a portfolio that demonstrates both theoretical understanding and practical skill.

The eight algorithms, in order of increasing complexity:

1. Bell State Preparation and Measurement
2. Deutsch-Jozsa Algorithm
3. Bernstein-Vazirani Algorithm
4. Quantum Fourier Transform
5. Grover's Search Algorithm
6. Simplified Shor's Factoring (N = 15)
7. Variational Quantum Eigensolver (VQE) for H₂
8. QAOA for MaxCut

For each algorithm, we follow the same structure: mathematical derivation → circuit diagram → Qiskit code → hardware execution → noise analysis.

The Gap Between Theory and Practice

Before we begin, let's set expectations. Every algorithm in this chapter works perfectly on a simulator. On real hardware, results will be imperfect due to noise. The gap between theoretical predictions and hardware results is not a failure — it is the central challenge of quantum computing. Understanding this gap quantitatively is more valuable than achieving perfect results on a simulator.

The noise budget for any quantum circuit has three main components:

1. Gate errors: Each gate operation has a finite fidelity. For a circuit with $d$ two-qubit gates at fidelity $F$, the overall circuit fidelity is approximately $F^d$.

2. Decoherence: Quantum states decay over time. For a circuit of duration $t$ on qubits with $T_1$ (energy relaxation) and $T_2$ (dephasing) times, the decoherence error is approximately $1 - e^{-t/T_{1,2}}$.

3. Readout errors: Measurement misclassification rates are typically 1-5% per qubit.

We will quantify these effects for each algorithm throughout this chapter.

33.2 Project 1: Bell State Preparation and Measurement

33.2.1 Mathematical Derivation

A Bell state is a maximally entangled two-qubit state. The four Bell states are:

$$|\Phi^+\rangle = \frac{1}{\sqrt{2}}(|00\rangle + |11\rangle), \quad |\Phi^-\rangle = \frac{1}{\sqrt{2}}(|00\rangle - |11\rangle),$$ $$|\Psi^+\rangle = \frac{1}{\sqrt{2}}(|01\rangle + |10\rangle), \quad |\Psi^-\rangle = \frac{1}{\sqrt{2}}(|01\rangle - |10\rangle).$$

The state $|\Phi^+\rangle$ is prepared by applying a Hadamard gate to qubit 0 followed by a CNOT with qubit 0 as control and qubit 1 as target:

$$|\Phi^+\rangle = \text{CNOT}_{0,1} \cdot (H_0 \otimes I_1) |00\rangle.$$

Proof:

$$|00\rangle \xrightarrow{H_0} \frac{1}{\sqrt{2}}(|0\rangle + |1\rangle) \otimes |0\rangle = \frac{1}{\sqrt{2}}(|00\rangle + |10\rangle)$$ $$\xrightarrow{\text{CNOT}_{0,1}} \frac{1}{\sqrt{2}}(|00\rangle + |11\rangle) = |\Phi^+\rangle.$$

Properties of Bell states:

Bell states are maximally entangled, meaning that measuring either qubit completely determines the other. For $|\Phi^+\rangle$:

• Measuring qubit 0 in the computational basis gives $|0\rangle$ or $|1\rangle$ with equal probability.
• If qubit 0 is measured as $|0\rangle$, qubit 1 is guaranteed to be $|0\rangle$.
• If qubit 0 is measured as $|1\rangle$, qubit 1 is guaranteed to be $|1\rangle$.

This perfect correlation holds regardless of the measurement basis. If both qubits are measured in the Hadamard basis ($|+\rangle, |-\rangle$), they remain perfectly correlated:

$$|\Phi^+\rangle = \frac{1}{\sqrt{2}}(|++\rangle + |--\rangle)$$

The reduced density matrix of either qubit is maximally mixed:

$$\rho_0 = \text{Tr}_1[|\Phi^+\rangle\langle\Phi^+|] = \frac{I}{2}$$

This means that knowing the state of one qubit tells you nothing about the other — until you measure it, at which point you know everything. This is the essence of entanglement.

Bell's inequality. The CHSH version of Bell's inequality provides a quantitative test:

$$S = E(a,b) + E(a,b') + E(a',b) - E(a',b') \leq 2$$

where $E(a,b)$ is the correlation between measurements in directions $a$ and $b$. Quantum mechanics predicts $S = 2\sqrt{2} \approx 2.828$ for optimally chosen angles, violating the classical bound. This violation has been experimentally confirmed to high precision, ruling out local hidden variable theories.

33.2.2 Circuit Diagram

ASCII Circuit: Bell State |Φ⁺⟩
=================================
q0: |0⟩ ──H──●── M (→ c0)
             │
q1: |0⟩ ─────X── M (→ c1)

Expected output: 00 with 50% probability, 11 with 50% probability.

33.2.3 Qiskit Implementation

from qiskit import QuantumCircuit, transpile
from qiskit_aer import AerSimulator
from qiskit.visualization import plot_histogram
import numpy as np

# ── Build Bell state circuit ──
def bell_state_circuit():
    qc = QuantumCircuit(2, 2)
    qc.h(0)
    qc.cx(0, 1)
    qc.measure([0, 1], [0, 1])
    return qc

qc_bell = bell_state_circuit()
print("Bell State Circuit:")
print(qc_bell.draw('text'))

# ── Simulate (ideal) ──
simulator = AerSimulator()
qc_bell_ideal = bell_state_circuit()
result_ideal = simulator.run(qc_bell_ideal, shots=8192).result()
counts_ideal = result_ideal.get_counts()

print("\nIdeal simulation results:")
for outcome, count in sorted(counts_ideal.items()):
    print(f"  |{outcome}⟩: {count} ({100*count/8192:.1f}%)")

# ── Run on real hardware ──
# Uncomment and configure with your IBM Quantum credentials:
#
# from qiskit_ibm_runtime import QiskitRuntimeService, Sampler
# service = QiskitRuntimeService(channel="ibm_quantum")
# backend = service.least_busy(simulator=False, operational=True)
# print(f"Selected backend: {backend.name}")
#
# qc_transpiled = transpile(qc_bell, backend=backend, optimization_level=3)
# job = backend.run(qc_transpiled, shots=8192)
# print(f"Job ID: {job.job_id}")
# result_hw = job.result()
# counts_hw = result_hw.get_counts()
#
# print("\nReal hardware results:")
# for outcome, count in sorted(counts_hw.items()):
#     print(f"  |{outcome}⟩: {count} ({100*count/8192:.1f}%)")

# ── Noise model simulation (approximates real hardware) ──
from qiskit_aer.noise import NoiseModel, depolarizing_error

noise_model = NoiseModel()
error_1q = depolarizing_error(0.001, 1)
error_2q = depolarizing_error(0.01, 2)
noise_model.add_all_qubit_quantum_error(error_1q, ['h', 'x'])
noise_model.add_all_qubit_quantum_error(error_2q, ['cx'])

result_noisy = simulator.run(
    qc_bell_ideal, noise_model=noise_model, shots=8192
).result()
counts_noisy = result_noisy.get_counts()

print("\nNoisy simulation results:")
for outcome, count in sorted(counts_noisy.items()):
    print(f"  |{outcome}⟩: {count} ({100*count/8192:.1f}%)")

# ── Compute fidelity ──
p_ideal_00 = counts_ideal.get('00', 0) / 8192
p_noisy_00 = counts_noisy.get('00', 0) / 8192
p_ideal_11 = counts_ideal.get('11', 0) / 8192
p_noisy_11 = counts_noisy.get('11', 0) / 8192

fidelity_estimate = (p_noisy_00 + p_noisy_11) / (p_ideal_00 + p_ideal_11)
print(f"\nEstimated Bell state fidelity: {fidelity_estimate:.4f}")

33.2.4 Noise Analysis

On real hardware, you will observe: - Leakage to $|01\rangle$ and $|10\rangle$: These outcomes should have zero probability ideally. Nonzero counts indicate gate errors (imperfect CNOT) or decoherence. - Asymmetry between $|00\rangle$ and $|11\rangle$: Amplitude damping (energy relaxation) preferentially drives $|11\rangle$ toward $|00\rangle$, creating an asymmetry. - Readout errors: Measurement misclassification flips some $|0\rangle$ outcomes to $|1\rangle$ and vice versa.

Quantitative noise budget for a Bell state circuit on typical IBM hardware:

Example 33.1: Computing Expected Fidelity

For a Bell state circuit on a specific backend with known calibration data:

def compute_expected_fidelity(t1_q0, t1_q1, t2_q0, t2_q1,
                              cx_error, h_error, meas_error_q0, meas_error_q1,
                              circuit_time):
    """
    Compute expected Bell state fidelity from hardware parameters.
    """
    # Gate errors
    gate_fidelity = (1 - h_error) * (1 - cx_error)

    # Decoherence
    t1_decay = np.exp(-circuit_time / min(t1_q0, t1_q1))
    t2_decay = np.exp(-circuit_time / min(t2_q0, t2_q1))
    decoherence_fidelity = 0.5 * (t1_decay + t2_decay)

    # Readout errors
    readout_fidelity = (1 - meas_error_q0) * (1 - meas_error_q1)

    # Combined (simplified model)
    total_fidelity = gate_fidelity * decoherence_fidelity * readout_fidelity

    return total_fidelity

# Typical IBM Brisbane parameters (2024)
fidelity = compute_expected_fidelity(
    t1_q0=300e-6, t1_q1=250e-6,      # T1 times in seconds
    t2_q0=200e-6, t2_q1=180e-6,      # T2 times in seconds
    cx_error=0.007,                     # 0.7% CNOT error
    h_error=0.0002,                     # 0.02% H error
    meas_error_q0=0.015,               # 1.5% readout error
    meas_error_q1=0.020,               # 2.0% readout error
    circuit_time=500e-9                  # 500 ns circuit time
)
print(f"Expected Bell state fidelity: {fidelity:.4f}")

Try It Yourself: Measure Bell State Fidelity on Real Hardware

Run the Bell state circuit on at least two different IBM Quantum backends. For each: 1. Record the backend name and its current calibration data (T1, T2, gate errors, readout errors). 2. Compute the expected fidelity from the noise budget model. 3. Measure the actual fidelity from the output counts. 4. Compare the expected and actual fidelities. Are they consistent? 5. Which backend gives higher fidelity? Does this correlate with the reported error rates?

33.3 Project 2: Deutsch-Jozsa Algorithm

33.3.1 Mathematical Derivation

The Deutsch-Jozsa algorithm determines whether a Boolean function $f: \{0,1\}^n \to \{0,1\}$ is constant (same output for all inputs) or balanced (outputs 0 for exactly half the inputs and 1 for the other half). Classically, this requires up to $2^{n-1} + 1$ queries in the worst case. Quantumly, it requires exactly one query.

The algorithm uses an oracle $U_f$ defined by:

$$U_f |x\rangle|y\rangle = |x\rangle|y \oplus f(x)\rangle.$$

The circuit is:

q0: |0⟩ ──H──●── ... ──●── H ── Measure
             │          │
q1: |0⟩ ──H──┼── ... ──┼── H ── Measure
             │          │
 ...         │   U_f    │
             │          │
qn: |0⟩ ──H──┼── ... ──┼── H ── Measure
             │          │
an: |1⟩ ──H──X── ... ──X── (ancilla, not measured)

Step-by-step derivation:

1. Initialize: $|0\rangle^{\otimes n}|1\rangle$

2. Apply Hadamard to all qubits:

$$\frac{1}{\sqrt{2^n}} \sum_{x=0}^{2^n-1} |x\rangle \otimes \frac{1}{\sqrt{2}}(|0\rangle - |1\rangle)$$

3. Apply the oracle $U_f$: The oracle maps $|x\rangle|y\rangle \to |x\rangle|y \oplus f(x)\rangle$. When the ancilla is in $|-\rangle = (|0\rangle - |1\rangle)/\sqrt{2}$, the result is a phase kickback:

$$\frac{1}{\sqrt{2^n}} \sum_{x=0}^{2^n-1} (-1)^{f(x)} |x\rangle \otimes |-\rangle$$

The ancilla qubit is unchanged! It serves only to apply the phase $(-1)^{f(x)}$ to each basis state $|x\rangle$.

4. Apply Hadamard to the first $n$ qubits: Using the identity $H^{\otimes n}|x\rangle = \frac{1}{\sqrt{2^n}} \sum_{z=0}^{2^n-1} (-1)^{x \cdot z} |z\rangle$ where $x \cdot z = \sum_i x_i z_i \pmod{2}$:

$$\frac{1}{2^n} \sum_{z=0}^{2^n-1} \left(\sum_{x=0}^{2^n-1} (-1)^{f(x) + x \cdot z}\right) |z\rangle \otimes |-\rangle$$

5. Measure: The probability of measuring $|0\rangle^{\otimes n}$ is:

$$P(|0\rangle^{\otimes n}) = \left| \frac{1}{2^n} \sum_{x} (-1)^{f(x)} \right|^2.$$

If $f$ is constant, $(-1)^{f(x)}$ is the same for all $x$, so $\sum_x (-1)^{f(x)} = \pm 2^n$, giving $P(|0\rangle^{\otimes n}) = 1$.

If $f$ is balanced, exactly half the terms are $+1$ and half are $-1$, so $\sum_x (-1)^{f(x)} = 0$, giving $P(|0\rangle^{\otimes n}) = 0$.

33.3.2 Qiskit Implementation (n = 3)

import numpy as np
from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator

def deutsch_jozsa_circuit(n, oracle_type='constant'):
    """
    Build the Deutsch-Jozsa circuit for n-bit input.

    Args:
        n: number of input qubits.
        oracle_type: 'constant' or 'balanced'.
    """
    qc = QuantumCircuit(n + 1, n)

    # Initialize ancilla to |1>
    qc.x(n)
    # Apply Hadamard to all qubits
    qc.h(range(n + 1))
    qc.barrier()

    # ── Oracle ──
    if oracle_type == 'constant':
        # Constant oracle: f(x) = 0 (do nothing) or f(x) = 1 (X on ancilla)
        # We choose f(x) = 0 for simplicity
        pass
    elif oracle_type == 'balanced':
        # Balanced oracle: f(x) = x_0 ⊕ x_1 ⊕ ... ⊕ x_{n-1}
        # Implement as CNOT from each input qubit to ancilla
        for i in range(n):
            qc.cx(i, n)

    qc.barrier()

    # Final Hadamard on input qubits
    qc.h(range(n))
    qc.measure(range(n), range(n))

    return qc

# ── Test ──
simulator = AerSimulator()

for oracle_type in ['constant', 'balanced']:
    qc = deutsch_jozsa_circuit(3, oracle_type)
    result = simulator.run(qc, shots=1024).result()
    counts = result.get_counts()
    print(f"\nDeutsch-Jozsa, {oracle_type} oracle:")
    print(f"  Result: {counts}")
    print(f"  All zeros? {'Yes — constant' if '000' in counts and len(counts) == 1 else 'No — balanced'}")

33.3.3 Noise Analysis

The Deutsch-Jozsa circuit for $n$ input qubits requires: - $(n+1)$ Hadamard gates (initialization) - $n$ CNOT gates (balanced oracle) - $n$ Hadamard gates (final layer) - Total: $(2n+1)$ single-qubit gates, $n$ two-qubit gates

For $n = 3$: 7 single-qubit gates, 3 two-qubit gates.

Expected fidelity on hardware with 0.02% single-qubit error and 0.7% two-qubit error:

$$F \approx (1 - 0.0002)^7 \times (1 - 0.007)^3 \times (1 - 0.015)^3 \approx 0.94$$

This means we expect ~94% success probability — the algorithm should work reliably on current hardware.

Common Misconception: "Deutsch-Jozsa proves quantum computers are faster."

Deutsch-Jozsa is an oracle separation: the quantum algorithm makes exactly one query, while the classical deterministic algorithm must make $2^{n-1}+1$ queries. But this is a worst-case classical bound, and the randomized classical algorithm can solve Deutsch-Jozsa with bounded error in only $O(1)$ queries (just query the oracle twice and check if both outputs are the same). The quantum advantage for Deutsch-Jozsa is exponential in the deterministic setting but only constant-factor in the randomized setting. It's an important conceptual demonstration, but not a practical speedup.

33.4 Project 3: Bernstein-Vazirani Algorithm

33.4.1 Mathematical Derivation

The Bernstein-Vazirani problem: given an oracle for $f(x) = s \cdot x \pmod{2}$ (the dot product of $x$ with a hidden string $s \in \{0,1\}^n$), find $s$. Classically, this requires $n$ queries (one per bit). Quantumly, it requires one query.

The circuit is identical to Deutsch-Jozsa, but the oracle encodes $s$:

$$U_f |x\rangle|y\rangle = |x\rangle|y \oplus (s \cdot x)\rangle.$$

After the oracle and final Hadamard gates, the measurement outcome is exactly $s$:

$$|\psi_{\text{final}}\rangle = |s\rangle \otimes |-\rangle.$$

Derivation:

After the oracle applies phase kickback $(-1)^{s \cdot x}$:

$$\frac{1}{\sqrt{2^n}} \sum_{x=0}^{2^n-1} (-1)^{s \cdot x} |x\rangle |-\rangle$$

After the final Hadamard:

$$H^{\otimes n} \left(\frac{1}{\sqrt{2^n}} \sum_{x=0}^{2^n-1} (-1)^{s \cdot x} |x\rangle\right) = |s\rangle$$

This follows from the identity $H^{\otimes n} \sum_x (-1)^{s \cdot x} |x\rangle = \sum_z \frac{1}{2^n} \sum_x (-1)^{s \cdot x + x \cdot z} |z\rangle$. The inner sum is:

$$\sum_{x=0}^{2^n-1} (-1)^{(s+z) \cdot x} = \begin{cases} 2^n & \text{if } s = z \\ 0 & \text{otherwise} \end{cases}$$

So the result is exactly $|s\rangle$ — the hidden string is revealed with certainty.

33.4.2 Qiskit Implementation

def bernstein_vazirani_circuit(n, hidden_string):
    """
    Build the Bernstein-Vazirani circuit.

    Args:
        n: number of input qubits.
        hidden_string: binary string of length n (the secret s).
    """
    qc = QuantumCircuit(n + 1, n)

    # Initialize ancilla to |1>
    qc.x(n)
    qc.h(range(n + 1))
    qc.barrier()

    # Oracle: for each bit i where s_i = 1, apply CNOT from qubit i to ancilla
    for i, bit in enumerate(hidden_string):
        if bit == '1':
            qc.cx(i, n)

    qc.barrier()
    qc.h(range(n))
    qc.measure(range(n), range(n))

    return qc

# ── Test ──
hidden = '10110'  # n = 5
qc_bv = bernstein_vazirani_circuit(5, hidden)
result = simulator.run(qc_bv, shots=1024).result()
counts = result.get_counts()
print(f"\nBernstein-Vazirani, hidden string s = {hidden}:")
print(f"  Measured: {list(counts.keys())[0]}")
print(f"  Correct:  {list(counts.keys())[0] == hidden}")

33.4.3 Noise Analysis

For a hidden string of weight $w$ (number of 1s), the Bernstein-Vazirani circuit uses $w$ CNOT gates plus $2(n+1)$ Hadamard gates. For $s = 10110$ (weight 3):

• 12 single-qubit gates, 3 two-qubit gates
• Expected fidelity: $\sim 96\%$
• This algorithm should work reliably on current hardware for $n \leq 8$.

Example 33.2: Scaling Analysis

How large can $n$ be before Bernstein-Vazirani fails on hardware? With single-qubit error rate $\epsilon_1 = 0.02\%$ and two-qubit error rate $\epsilon_2 = 0.7\%$:

• Total single-qubit gates: $2n + 2$
• Total two-qubit gates: $\text{weight}(s) \leq n$
• Expected fidelity: $(1 - \epsilon_1)^{2n+2} \cdot (1 - \epsilon_2)^{n}$
• Setting $F = 0.5$ (failure threshold): solve for $n$

$$n \approx \frac{\ln 0.5}{\ln(1-\epsilon_1) \cdot 2 + \ln(1-\epsilon_2)} \approx 70$$

So Bernstein-Vazirani should work reliably up to ~70 bits on current hardware — well beyond the classical $n$-query bound. This demonstrates the practical advantage of the quantum approach.

33.5 Project 4: Quantum Fourier Transform

33.5.1 Mathematical Derivation

The Quantum Fourier Transform (QFT) on $n$ qubits maps a computational basis state $|j\rangle$ to:

$$\text{QFT}|j\rangle = \frac{1}{\sqrt{2^n}} \sum_{k=0}^{2^n-1} e^{2\pi i j k / 2^n} |k\rangle.$$

For an arbitrary state $|\psi\rangle = \sum_j x_j |j\rangle$, the QFT produces:

$$\text{QFT}|\psi\rangle = \sum_k \tilde{x}_k |k\rangle, \quad \tilde{x}_k = \frac{1}{\sqrt{2^n}} \sum_{j=0}^{2^n-1} x_j e^{2\pi i j k / 2^n}.$$

This is the discrete Fourier transform of the amplitudes — the quantum analog of the classical FFT, but operating on amplitudes rather than data.

The key difference from the classical FFT: the QFT operates on the amplitudes of a quantum state, not on classical data. You cannot read out all $2^n$ amplitudes — measuring gives you one basis state. This is why the QFT by itself doesn't provide exponential speedup for classical Fourier transform; it's useful as a subroutine in algorithms like Shor's, where the periodic structure of the measurement outcome is what matters.

The QFT circuit uses controlled phase rotations $R_k$:

$$R_k = \begin{pmatrix} 1 & 0 \\ 0 & e^{2\pi i / 2^k} \end{pmatrix}.$$

The circuit for $n$ qubits requires $\frac{n(n-1)}{2}$ controlled-$R_k$ gates plus $n$ Hadamard gates, for a total depth of $O(n^2)$. With parallelization, this can be reduced to $O(n)$ depth.

33.5.2 Circuit Diagram (n = 3)

ASCII Circuit: QFT on 3 Qubits
================================
q0: ──H──R2──R3─────────────────────●──────●──
         │   │                       │      │
q1: ─────●───┼──H──R2────────────────●──────┼──
             │   │                  │      │
q2: ─────────●───●──H───────────────┼──────●──
                                    │
Swap q0 and q2 for standard ordering.

33.5.3 Qiskit Implementation

import numpy as np
from qiskit import QuantumCircuit
from qiskit.circuit.library import QFT
from qiskit_aer import AerSimulator
from qiskit.quantum_info import Statevector

def qft_demo(n_qubits=3):
    """Demonstrate the QFT on a specific input state."""
    # ── Build QFT circuit ──
    qc = QuantumCircuit(n_qubits)

    # Prepare input state |j> = |5> = |101> for n=3
    qc.x(0)  # Set qubit 0 to |1>
    qc.x(2)  # Set qubit 2 to |1>
    # State is now |101> = |5>

    qc.barrier()

    # Apply QFT
    qft = QFT(num_qubits=n_qubits, do_swaps=True)
    qc.append(qft, range(n_qubits))

    # ── Get statevector ──
    state = Statevector.from_instruction(qc)
    print(f"\nQFT of |5> (n={n_qubits}):")
    for i, amp in enumerate(state.data):
        if abs(amp) > 1e-10:
            print(f"  |{i:0{n_qubits}b}⟩: {amp:.4f}")

    # ── Verify: QFT|j> = (1/√N) Σ_k exp(2πi j k / N) |k> ──
    N = 2**n_qubits
    j = 5
    print("\nTheoretical amplitudes:")
    for k in range(N):
        amp_theory = np.exp(2j * np.pi * j * k / N) / np.sqrt(N)
        if abs(amp_theory) > 1e-10:
            print(f"  |{k:0{n_qubits}b}⟩: {amp_theory:.4f}")

    return qc

qc_qft = qft_demo(3)
print("\nQFT Circuit:")
print(qc_qft.draw('text'))

# ── Inverse QFT ──
def qft_inverse_demo():
    """Demonstrate QFT followed by inverse QFT returns original state."""
    qc = QuantumCircuit(3)
    qc.x(0)
    qc.x(2)  # |101>

    qft = QFT(num_qubits=3, do_swaps=True)
    iqft = QFT(num_qubits=3, do_swaps=True, inverse=True)

    qc.append(qft, range(3))
    qc.append(iqft, range(3))

    state = Statevector.from_instruction(qc)
    print("\nQFT⁻¹(QFT(|5>)):")
    for i, amp in enumerate(state.data):
        if abs(amp) > 1e-10:
            print(f"  |{i:03b}⟩: {amp:.4f}")

qft_inverse_demo()

33.5.4 Noise Analysis

The QFT on $n$ qubits uses $\frac{n(n-1)}{2}$ controlled-$R_k$ gates. For $n = 3$: 3 controlled rotations + 3 Hadamard gates + 2 SWAP gates.

On hardware, controlled rotations are decomposed into CNOT gates and single-qubit rotations. Each controlled-$R_k$ requires 2 CNOTs. Total CNOT count: $n(n-1) + n$ (from SWAPs) $\approx n^2$.

For $n = 3$: ~9 CNOT gates. Expected fidelity: $\sim 93\%$. For $n = 5$: ~25 CNOT gates. Expected fidelity: $\sim 83\%$. For $n = 10$: ~100 CNOT gates. Expected fidelity: $\sim 50\%$.

The QFT becomes unreliable on NISQ hardware for $n > 10$ without error mitigation.

Try It Yourself: QFT Fidelity Scaling

Implement the QFT for $n = 2, 3, 4, 5, 6$ and measure the fidelity of the inverse-QFT test (QFT followed by inverse QFT should return the original state) on a noisy simulator. Plot the fidelity as a function of $n$. At what $n$ does the fidelity drop below 90%? Below 50%? How does this compare with the theoretical prediction $F \approx (1 - \epsilon_{CX})^{n^2}$?

33.6 Project 5: Grover's Search Algorithm

33.6.1 Mathematical Derivation

Grover's algorithm searches an unstructured database of $N = 2^n$ items for a marked element $m$ using $O(\sqrt{N})$ queries. The algorithm iterates a Grover operator $G$:

$$G = H^{\otimes n} (2|0\rangle\langle 0|^{\otimes n} - I) H^{\otimes n} \cdot O,$$

where $O$ is the oracle that flips the phase of the marked state:

$$O|x\rangle = (-1)^{f(x)}|x\rangle, \quad f(x) = \begin{cases} 1 & x = m \\ 0 & x \neq m \end{cases}.$$

Derivation of the optimal number of iterations:

The state after $k$ Grover iterations is approximately:

$$|\psi_k\rangle \approx \sin\left(\frac{(2k+1)\theta}{2}\right) |m\rangle + \cos\left(\frac{(2k+1)\theta}{2}\right) |m^\perp\rangle$$

where $\sin\theta = \sqrt{M/N}$ for $M$ marked items (here $M = 1$). The probability of measuring the marked state is $\sin^2((2k+1)\theta/2)$, which is maximized when $(2k+1)\theta/2 = \pi/2$, giving:

$$k_{\text{opt}} \approx \frac{\pi}{4\theta} \approx \frac{\pi}{4}\sqrt{\frac{N}{M}}$$

For $M = 1$: $k_{\text{opt}} \approx \lfloor \pi\sqrt{N}/4 \rfloor$.

33.6.2 Qiskit Implementation (n = 3, search for |101⟩)

import numpy as np
from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator
from qiskit.circuit.library import GroverOperator, MCMT, ZGate

def grover_search(n_qubits, marked_state):
    """
    Grover's search for a specific marked state.

    Args:
        n_qubits: number of qubits (database size = 2^n_qubits).
        marked_state: binary string of the marked element.
    """
    qc = QuantumCircuit(n_qubits, n_qubits)

    # ── Initialize superposition ──
    qc.h(range(n_qubits))
    qc.barrier()

    # ── Number of iterations ──
    N = 2**n_qubits
    iterations = int(np.floor(np.pi / 4 * np.sqrt(N)))

    # ── Build oracle for the marked state ──
    def oracle(marked):
        """Phase oracle that flips the phase of |marked>."""
        oc = QuantumCircuit(n_qubits)
        marked_bits = [int(b) for b in marked]
        for i, b in enumerate(marked_bits):
            if b == 0:
                oc.x(i)
        oc.h(n_qubits - 1)
        oc.mcx(list(range(n_qubits - 1)), n_qubits - 1)
        oc.h(n_qubits - 1)
        for i, b in enumerate(marked_bits):
            if b == 0:
                oc.x(i)
        return oc

    # ── Grover iterations ──
    for _ in range(iterations):
        # Oracle
        qc.append(oracle(marked_state), range(n_qubits))
        qc.barrier()
        # Diffusion operator
        qc.h(range(n_qubits))
        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))
        qc.h(range(n_qubits))
        qc.barrier()

    qc.measure(range(n_qubits), range(n_qubits))
    return qc, iterations

# ── Run ──
marked = '101'
qc_grover, iters = grover_search(3, marked)
print(f"Grover's Search for |{marked}⟩, {iters} iterations:")
print(qc_grover.draw('text'))

simulator = AerSimulator()
result = simulator.run(qc_grover, shots=8192).result()
counts = result.get_counts()
print(f"\nResults:")
for outcome, count in sorted(counts.items(), key=lambda x: -x[1]):
    marker = " ← MARKED" if outcome == marked else ""
    print(f"  |{outcome}⟩: {count:>4} ({100*count/8192:.1f}%){marker}")

33.6.3 Noise Analysis

For $n = 3$ qubits with 1 marked item, Grover requires 2 iterations. Each iteration uses: - Oracle: 2 X gates + 1 multi-controlled Z (decomposed into 6 CNOTs) + 2 X gates - Diffusion: same gate count as oracle - Total per iteration: ~16 CNOTs + ~24 single-qubit gates - Total for 2 iterations: ~32 CNOTs + ~48 single-qubit gates

Expected fidelity: $(1 - 0.007)^{32} \times (1 - 0.0002)^{48} \approx 0.80$

This means about 80% of shots will give the correct answer, with 20% distributed among wrong answers. On current hardware, Grover's search works for small databases ($n \leq 4$) but becomes unreliable for larger ones.

Example 33.3: Grover's Search for Multiple Marked Items

When there are $M$ marked items out of $N$, the optimal number of Grover iterations changes:

$$k_{\text{opt}} = \left\lfloor \frac{\pi}{4} \sqrt{\frac{N}{M}} \right\rfloor$$

The success probability after $k$ iterations is $\sin^2((2k+1)\arcsin\sqrt{M/N})$.

# Compute success probability as a function of iterations
def grover_success_prob(N, M, k):
    """Success probability after k Grover iterations."""
    theta = np.arcsin(np.sqrt(M / N))
    return np.sin((2 * k + 1) * theta) ** 2

# For N=8, M=1 (1 marked out of 8)
N, M = 8, 1
for k in range(5):
    p = grover_success_prob(N, M, k)
    print(f"  k={k}: P(success) = {p:.4f}")

# Optimal: k = floor(pi/4 * sqrt(8)) = 2
# Success probability: ~0.9453

33.7 Project 6: Simplified Shor's Factoring (N = 15)

33.7.1 Mathematical Derivation

Shor's algorithm factors $N = 15$ by finding the period $r$ of $f(x) = a^x \bmod 15$. We choose $a = 2$ (or $a = 7, 8, 11, 13$ — all coprime to 15). The function values are:

$$f(0) = 1, \; f(1) = 2, \; f(2) = 4, \; f(3) = 8, \; f(4) = 1, \; \ldots$$

The period is $r = 4$. Then $\gcd(a^{r/2} \pm 1, N) = \gcd(2^2 \pm 1, 15) = \gcd(3, 15) = 3$ and $\gcd(5, 15) = 5$, yielding the factors.

Detailed step-by-step for $a = 7$:

1. Choose $a = 7$ (coprime to 15, since $\gcd(7, 15) = 1$).
2. Compute $f(x) = 7^x \bmod 15$: - $7^1 \bmod 15 = 7$ - $7^2 \bmod 15 = 49 \bmod 15 = 4$ - $7^3 \bmod 15 = 343 \bmod 15 = 13$ - $7^4 \bmod 15 = 2401 \bmod 15 = 1$
3. Period $r = 4$.
4. Check: $r$ is even, so we can proceed.
5. Compute: $a^{r/2} = 7^2 = 49$.
6. Factors: $\gcd(49 - 1, 15) = \gcd(48, 15) = 3$ and $\gcd(49 + 1, 15) = \gcd(50, 15) = 5$.
7. Result: $15 = 3 \times 5$. ✓

Why does this work?

If $a^r \equiv 1 \pmod{N}$, then $a^r - 1 \equiv 0 \pmod{N}$, which means $N | (a^{r/2} - 1)(a^{r/2} + 1)$. If $a^{r/2} \not\equiv \pm 1 \pmod{N}$, then both $a^{r/2} - 1$ and $a^{r/2} + 1$ have non-trivial common factors with $N$.

33.7.2 Qiskit Implementation

import numpy as np
from math import gcd
from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator
from qiskit.circuit.library import QFT

def shor_15(a=7):
    """
    Simplified Shor's algorithm for N=15.
    Uses a pre-optimized circuit for modular exponentiation.

    Args:
        a: base coprime to 15 (2, 7, 8, 11, 13).
    """
    # ── Circuit: 8 qubits for period estimation, 4 for modular exponentiation ──
    n_count = 8
    n_mod = 4
    total = n_count + n_mod

    qc = QuantumCircuit(total, n_count)

    # Initialize modular register to |1>
    qc.x(n_count)  # Set the first qubit of modular register to |1>

    # Superposition on counting register
    qc.h(range(n_count))
    qc.barrier()

    # ── Modular exponentiation: controlled a^(2^j) mod 15 ──
    # For a=7: 7^1=7, 7^2=4, 7^4=1, 7^8=1, ...
    # We implement controlled multiplications using SWAP gates
    # (This is the simplified circuit for N=15, a=7)

    # Controlled-U^(2^0): multiply by 7 mod 15
    qc.cswap(n_count - 1, n_count + 1, n_count + 2)
    qc.cswap(n_count - 1, n_count + 2, n_count + 3)
    qc.cswap(n_count - 1, n_count + 0, n_count + 3)

    # Controlled-U^(2^1): multiply by 4 mod 15
    qc.cswap(n_count - 2, n_count + 0, n_count + 2)
    qc.cswap(n_count - 2, n_count + 1, n_count + 3)

    # Controlled-U^(2^2): multiply by 1 mod 15 (identity, skip)
    # Controlled-U^(2^3): multiply by 1 mod 15 (identity, skip)

    qc.barrier()

    # ── Inverse QFT on counting register ──
    iqft = QFT(num_qubits=n_count, inverse=True, do_swaps=True)
    qc.append(iqft, range(n_count))

    qc.measure(range(n_count), range(n_count))
    return qc

# ── Run ──
qc_shor = shor_15(a=7)
print("Shor's Algorithm for N=15, a=7:")
print(qc_shor.draw('text'))

simulator = AerSimulator()
result = simulator.run(qc_shor, shots=8192).result()
counts = result.get_counts()

print("\nMeasurement results (top 8):")
for outcome, count in sorted(counts.items(), key=lambda x: -x[1])[:8]:
    phase = int(outcome, 2) / 256  # 2^8 = 256
    print(f"  |{outcome}⟩: {count:>4} shots, phase = {phase:.4f}")

# ── Classical post-processing ──
def find_period_from_measurements(counts, N=15, a=7):
    """Extract period from measurement results using continued fractions."""
    from fractions import Fraction

    total_shots = sum(counts.values())
    phases = {}

    for bitstring, count in counts.items():
        measured = int(bitstring, 2)
        phase = measured / 256
        # Find closest fraction
        frac = Fraction(phase).limit_denominator(15)
        r_candidate = frac.denominator

        if r_candidate > 0 and pow(a, r_candidate, N) == 1:
            phases[r_candidate] = phases.get(r_candidate, 0) + count

    if not phases:
        return None

    best_r = max(phases, key=phases.get)
    return best_r

r = find_period_from_measurements(counts)
print(f"\nEstimated period r = {r}")

if r is not None and r % 2 == 0:
    factor1 = gcd(7**(r//2) - 1, 15)
    factor2 = gcd(7**(r//2) + 1, 15)
    print(f"Factors of 15: {factor1} × {factor2}")

33.7.3 Noise Analysis

Shor's algorithm for N=15 uses 12 qubits and approximately 50 CNOT gates (after decomposition). This is at the limit of what current NISQ hardware can handle.

• Expected fidelity: $(1 - 0.007)^{50} \approx 0.70$ on hardware with 0.7% CNOT error
• On real hardware: The circuit typically works with ~30-50% success probability on the best IBM backends, requiring multiple shots to extract the correct period.

The circuit depth is too high for reliable execution on current hardware. This illustrates a key theme: algorithms that require deep circuits (like Shor's) are not yet feasible on NISQ devices, while shallow circuits (like Bell states) work well.

Common Misconception: "Shor's algorithm factoring 15 proves quantum computers can break RSA."

Factoring 15 is a proof of principle, not a practical demonstration. The circuit for N=15 is heavily optimized and specific to that number. Factoring RSA-2048 would require ~20 million physical qubits and hours of computation time. The N=15 demonstration shows that the quantum period-finding subroutine works, but it doesn't demonstrate cryptographically relevant factoring.

33.8 Project 7: VQE for H₂ Molecule

33.8.1 Mathematical Derivation

The Variational Quantum Eigensolver (VQE) finds the ground state energy of a Hamiltonian $H$ by minimizing:

$$E(\boldsymbol{\theta}) = \langle \psi(\boldsymbol{\theta}) | H | \psi(\boldsymbol{\theta}) \rangle,$$

where $|\psi(\boldsymbol{\theta})\rangle = U(\boldsymbol{\theta})|0\rangle$ is a parameterized trial state. The Hamiltonian is expressed as a sum of Pauli strings:

$$H = \sum_i c_i P_i, \quad P_i \in \{I, X, Y, Z\}^{\otimes n}.$$

For H₂ in the minimal STO-3G basis, the qubit Hamiltonian (after Jordan-Wigner or Bravyi-Kitaev transformation) has 4 qubits and 15 terms.

The variational principle guarantees that $E(\boldsymbol{\theta}) \geq E_0$ (the true ground state energy) for any choice of $\boldsymbol{\theta}$. The VQE algorithm:

1. Prepare the ansatz state $|\psi(\boldsymbol{\theta})\rangle$ on a quantum computer.
2. Measure each Pauli term $\langle P_i \rangle$ separately.
3. Compute $E(\boldsymbol{\theta}) = \sum_i c_i \langle P_i \rangle$.
4. Use a classical optimizer to update $\boldsymbol{\theta}$.
5. Repeat until convergence.

The choice of ansatz is critical: - Hardware-efficient ansatz: Alternating layers of single-qubit rotations and entangling gates. Expressive but prone to barren plateaus. - UCCSD (Unitary Coupled Cluster Singles and Doubles): Chemically motivated ansatz that respects particle number and spin symmetry. Deeper circuits but better convergence. - Adaptive ansatz (ADAPT-VQE): Grow the circuit iteratively by adding operators that reduce the energy the most.

33.8.2 Qiskit Implementation

import numpy as np
from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator
from qiskit.quantum_info import SparsePauliOp
from qiskit.circuit.library import EfficientSU2
from scipy.optimize import minimize

# ── H₂ Hamiltonian (STO-3G, bond distance 0.735 Å) ──
def h2_hamiltonian():
    """H₂ Hamiltonian in the minimal basis, 4 qubits."""
    coeffs = [
        -0.810547,   # II
        0.172183,    # Z0
        0.172183,    # Z1
        -0.225753,   # Z2
        -0.225753,   # Z3
        0.120912,    # Z0 Z1
        0.168927,    # Z0 Z2
        0.166145,    # Z0 Z3
        0.166145,    # Z1 Z2
        0.168927,    # Z1 Z3
        0.174643,    # Z2 Z3
        0.045232,    # Y0 Y1 X2 X3
        0.045232,    # X0 X1 Y2 Y3
        0.045232,    # Y0 X1 X2 Y3
        0.045232,    # X0 Y1 Y2 X3
    ]
    paulis = [
        'IIII', 'ZIII', 'IZII', 'IIZI', 'IIIZ',
        'ZZII', 'ZIZI', 'ZIIZ', 'IZZI', 'IZIZ', 'IIZZ',
        'YYXX', 'XXYY', 'YXXY', 'XYYX',
    ]
    return SparsePauliOp(paulis, coeffs)

H = h2_hamiltonian()
print(f"H₂ Hamiltonian: {len(H)} terms")

# ── Ansatz circuit ──
def build_ansatz(params, n_qubits=4, reps=2):
    """Build the parameterized ansatz circuit."""
    ansatz = EfficientSU2(n_qubits, su2_gates=['ry', 'rz'],
                          entanglement='full', reps=reps)
    qc = QuantumCircuit(n_qubits)
    qc.compose(ansatz.assign_parameters(params), inplace=True)
    return qc

# ── Energy evaluation ──
def compute_energy(params, H, n_qubits=4, reps=2):
    """Compute <ψ(θ)|H|ψ(θ)> using the AerSimulator."""
    qc = build_ansatz(params, n_qubits, reps)
    simulator = AerSimulator(method='statevector')
    result = simulator.run(qc).result()
    state = result.get_statevector()
    energy = state.expectation_value(H).real
    return energy

# ── VQE optimization ──
n_qubits = 4
reps = 2
ansatz_blank = EfficientSU2(n_qubits, su2_gates=['ry', 'rz'],
                             entanglement='full', reps=reps)
n_params = ansatz_blank.num_parameters

print(f"Number of parameters: {n_params}")

# Initial parameters
np.random.seed(42)
initial_params = np.random.uniform(0, 2*np.pi, n_params)

# Run optimization
print("Running VQE optimization...")
result = minimize(
    compute_energy,
    initial_params,
    args=(H, n_qubits, reps),
    method='COBYLA',
    options={'maxiter': 200, 'disp': True}
)

print(f"\nVQE Results:")
print(f"  Optimized energy: {result.fun:.6f} Hartree")
print(f"  Exact ground state energy (H₂ at 0.735 Å): -1.857275 Hartree")
print(f"  Error: {abs(result.fun - (-1.857275)):.6f} Hartree")

33.8.3 Noise Analysis and Error Mitigation

VQE is relatively resilient to noise because the variational principle bounds the energy: noise can only increase the estimated energy (for a reasonable ansatz), and the optimizer can partially compensate by finding parameters that minimize the noisy energy.

Zero-noise extrapolation (ZNE) is a powerful error mitigation technique:

1. Run the circuit at the original noise level $\lambda = 1$.
2. Run the circuit at amplified noise levels $\lambda = 2, 3$ (by stretching gate durations or inserting identity pairs).
3. Extrapolate the results to $\lambda = 0$ (zero noise).

def zero_noise_extrapolation(energy_at_lambda_1, energy_at_lambda_2, energy_at_lambda_3):
    """
    Richardson extrapolation from three noise levels.
    lambda=1: original circuit
    lambda=2: each gate repeated twice (or pulse stretched by 2x)
    lambda=3: each gate repeated three times
    """
    # Linear extrapolation
    # E(0) = 3*E(1) - 3*E(2) + E(3) (Richardson extrapolation)
    return 3 * energy_at_lambda_1 - 3 * energy_at_lambda_2 + energy_at_lambda_3

Example 33.4: VQE with Noise

Running VQE on a realistic noise model for an IBM Quantum processor:

This demonstrates that with error mitigation, VQE can achieve chemical accuracy for small molecules even on noisy hardware.

33.9 Project 8: QAOA for MaxCut

33.9.1 Mathematical Derivation

The Quantum Approximate Optimization Algorithm (QAOA) solves combinatorial optimization problems. For MaxCut on a graph $G = (V, E)$, the goal is to partition vertices into two sets to maximize the number of edges between them.

The cost Hamiltonian is:

$$H_C = \frac{1}{2} \sum_{(i,j) \in E} (I - Z_i Z_j).$$

The QAOA circuit alternates between the cost unitary $U_C(\gamma) = e^{-i\gamma H_C}$ and the mixer unitary $U_B(\beta) = e^{-i\beta \sum_i X_i}$:

$$|\psi(\boldsymbol{\gamma}, \boldsymbol{\beta})\rangle = U_B(\beta_p) U_C(\gamma_p) \cdots U_B(\beta_1) U_C(\gamma_1) |+\rangle^{\otimes n}.$$

Performance guarantee: For $p = 1$ QAOA on MaxCut for 3-regular graphs, Farhi et al. proved that the expected cut value is at least 0.692 × optimal. This is better than random guessing (0.5 × optimal) but worse than the best classical algorithm (Goemans-Williamson: 0.878 × optimal).

The QAOA landscape (cost as a function of $\gamma$ and $\beta$) is generally non-convex with many local minima. For $p = 1$, the landscape can be visualized as a 2D surface, and the global minimum can be found by grid search. For $p \geq 2$, the landscape has exponentially many local minima, making optimization challenging.

33.9.2 Qiskit Implementation

import numpy as np
from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator
from scipy.optimize import minimize

def maxcut_qaoa(adjacency_matrix, p=2):
    """
    QAOA for MaxCut on a graph.

    Args:
        adjacency_matrix: n x n symmetric matrix (0/1 for unweighted).
        p: number of QAOA layers.
    """
    n = len(adjacency_matrix)

    def build_qaoa_circuit(gamma, beta):
        qc = QuantumCircuit(n, n)

        # Initial state: |+>^⊗n
        qc.h(range(n))

        for layer in range(p):
            # Cost unitary: e^{-i γ H_C}
            for i in range(n):
                for j in range(i + 1, n):
                    if adjacency_matrix[i, j] != 0:
                        # e^{-i γ (I - Z_i Z_j)/2} = e^{-iγ/2} * e^{iγ Z_i Z_j / 2}
                        # Implement as: CNOT(i,j), Rz(γ) on j, CNOT(i,j)
                        qc.cx(i, j)
                        qc.rz(2 * gamma[layer], j)
                        qc.cx(i, j)

            # Mixer unitary: e^{-i β Σ X_i}
            for i in range(n):
                qc.rx(2 * beta[layer], i)

        qc.measure(range(n), range(n))
        return qc

    def compute_expectation(params):
        gamma = params[:p]
        beta = params[p:]

        qc = build_qaoa_circuit(gamma, beta)
        simulator = AerSimulator()
        result = simulator.run(qc, shots=4096).result()
        counts = result.get_counts()

        # Compute expected cut value
        avg_cut = 0
        total_shots = sum(counts.values())
        for bitstring, count in counts.items():
            # Each bit is a partition assignment
            cut_value = 0
            for i in range(n):
                for j in range(i + 1, n):
                    if adjacency_matrix[i, j] != 0:
                        if bitstring[i] != bitstring[j]:
                            cut_value += adjacency_matrix[i, j]
            avg_cut += cut_value * count / total_shots

        return -avg_cut  # Negative because we minimize

    # Initial parameters
    np.random.seed(42)
    initial_params = np.random.uniform(0, np.pi, 2 * p)

    result = minimize(
        compute_expectation,
        initial_params,
        method='COBYLA',
        options={'maxiter': 100}
    )

    return result

# ── Test on a 4-node graph (square) ──
adj = np.array([
    [0, 1, 0, 1],
    [1, 0, 1, 0],
    [0, 1, 0, 1],
    [1, 0, 1, 0],
])

print("QAOA for MaxCut on a 4-node cycle graph")
result = maxcut_qaoa(adj, p=2)
print(f"  Optimal parameters: {result.x}")
print(f"  Max cut value: {-result.fun:.2f}")
print(f"  Optimal cut: 4 (all edges cut)")

33.9.3 Noise Analysis

QAOA is one of the most NISQ-friendly algorithms because: 1. The circuit depth is shallow (proportional to $p \times |E|$, where $|E|$ is the number of edges). 2. The variational nature allows the optimizer to partially compensate for noise. 3. The output is a classical bitstring, so readout errors are the main concern.

For the 4-node cycle graph with $p = 2$: - Circuit depth: ~16 CNOT gates + ~16 RX gates + measurements - Expected fidelity on hardware: ~85-90% - The algorithm typically produces the correct answer in ~60-70% of shots, which is sufficient for most optimization tasks.

Try It Yourself: QAOA Parameter Landscape

For a 3-regular graph with 6 vertices, plot the QAOA cost function landscape for $p = 1$ as a function of $\gamma$ and $\beta$. Use a grid of $50 \times 50$ points. Where are the local minima? How does the landscape change for $p = 2$ (4 parameters)? Can you identify the global minimum?

33.10 Running on Real Hardware

33.10.1 IBM Quantum Setup

To run any of the above algorithms on real IBM Quantum hardware:

from qiskit_ibm_runtime import QiskitRuntimeService, Sampler, Session
from qiskit import transpile

# 1. Authenticate
service = QiskitRuntimeService(
    channel="ibm_quantum",
    token="YOUR_IBM_QUANTUM_TOKEN"  # Get from https://quantum-computing.ibm.com
)

# 2. Select a backend
backend = service.least_busy(
    simulator=False,
    operational=True,
    min_num_qubits=2
)
print(f"Selected backend: {backend.name}")
print(f"  Qubits: {backend.configuration().num_qubits}")

# 3. Transpile for the target backend
qc = bell_state_circuit()  # Or any other circuit
qc_transpiled = transpile(qc, backend=backend, optimization_level=3)

# 4. Run using the Sampler primitive
with Session(service=service, backend=backend) as session:
    sampler = Sampler(session=session)
    job = sampler.run([qc_transpiled], shots=4096)
    print(f"Job ID: {job.job_id()}")

# 5. Retrieve results
result = job.result()
quasi_dists = result.quasi_dists[0]
print(f"Quasi-distribution: {quasi_dists}")

33.10.2 Interpreting Hardware Results

When comparing ideal and hardware results, compute:

1. Hellinger fidelity: $F_H = \left( \sum_x \sqrt{p_{\text{ideal}}(x) p_{\text{hardware}}(x)} \right)^2$
2. Total variation distance: $D_{TV} = \frac{1}{2} \sum_x |p_{\text{ideal}}(x) - p_{\text{hardware}}(x)|$
3. Success probability: Fraction of shots yielding the correct answer.

def compare_results(counts_ideal, counts_hardware):
    """Compare ideal and hardware measurement results."""
    total_ideal = sum(counts_ideal.values())
    total_hw = sum(counts_hardware.values())

    all_outcomes = set(counts_ideal.keys()) | set(counts_hardware.keys())

    p_ideal = {k: counts_ideal.get(k, 0) / total_ideal for k in all_outcomes}
    p_hw = {k: counts_hardware.get(k, 0) / total_hw for k in all_outcomes}

    # Hellinger fidelity
    hellinger = sum(np.sqrt(p_ideal[k] * p_hw[k]) for k in all_outcomes) ** 2

    # Total variation distance
    tv = 0.5 * sum(abs(p_ideal[k] - p_hw[k]) for k in all_outcomes)

    print(f"Hellinger fidelity: {hellinger:.4f}")
    print(f"Total variation distance: {tv:.4f}")

    return hellinger, tv

33.10.3 Error Mitigation Techniques

Beyond zero-noise extrapolation, several error mitigation techniques can improve hardware results:

Measurement error mitigation: Calibrate the readout error matrix and apply its inverse to the measured counts.

from qiskit.result import LocalReadoutMitigator

# Generate calibration circuits
cal_circuits = []
for state in ['00', '01', '10', '11']:
    qc = QuantumCircuit(2, 2)
    for i, bit in enumerate(state):
        if bit == '1':
            qc.x(i)
    qc.measure([0, 1], [0, 1])
    cal_circuits.append(qc)

# Run calibration circuits and compute mitigator
# (simplified; in practice, use Qiskit's complete measurement mitigation)

Dynamical decoupling: Insert identity-equivalent pulse sequences during idle qubit periods to reduce decoherence.

from qiskit import transpile

# Transpile with dynamical decoupling
qc_dd = transpile(qc, backend=backend,
                  optimization_level=3,
                  scheduling_method='asap',
                  instruction_durations=backend.instruction_durations)

33.11 Portfolio Summary

Key observations from the portfolio:

1. Shallow circuits work well. Algorithms with fewer than ~20 CNOT gates produce reliable results on current hardware.
2. Deep circuits struggle. Shor's algorithm for N=15 is at the edge of feasibility.
3. Variational algorithms are resilient. VQE and QAOA can tolerate noise because the optimizer compensates.
4. Error mitigation helps. Zero-noise extrapolation and readout error mitigation can improve results by 5-10×.
5. The gap between theory and practice is measurable and quantifiable. This is the central lesson of the capstone.

33.12 Advanced Topic: Noise-Aware Circuit Compilation

When running quantum algorithms on real hardware, the choice of qubit mapping and gate decomposition significantly affects performance. This section covers the key techniques for optimizing circuits for noisy hardware.

Qubit Mapping and Routing

On a real quantum processor, qubits are arranged in a specific connectivity topology (e.g., heavy-hex for IBM processors). A quantum circuit that assumes all-to-all connectivity must be mapped to this topology, requiring SWAP gates to move quantum information between non-adjacent qubits.

Each SWAP gate costs 3 CNOT gates, adding significant error. The transpiler must find a mapping that minimizes the total number of SWAPs:

from qiskit import transpile
from qiskit_ibm_runtime import QiskitRuntimeService

# Connect to IBM Quantum
service = QiskitRuntimeService(channel="ibm_quantum")
backend = service.least_busy(simulator=False, operational=True)

# Transpile the Bell state circuit for this backend
qc_bell_mapped = transpile(
    bell_state_circuit(),
    backend=backend,
    optimization_level=3,  # Maximum optimization
    routing_method='sabre',  # Best routing algorithm
    layout_method='sabre'   # Best layout algorithm
)

print(f"Original circuit depth: {bell_state_circuit().depth()}")
print(f"Transpiled circuit depth: {qc_bell_mapped.depth()}")
print(f"Original CNOT count: {bell_state_circuit().count_ops().get('cx', 0)}")
print(f"Transpiled CNOT count: {qc_bell_mapped.count_ops().get('cx', 0)}")

Gate Decomposition

The hardware-native gate set differs from the abstract gate set. For example, IBM processors natively support:

• Single-qubit: $\sqrt{X}$ (SX), $R_z(\theta)$, $X$
• Two-qubit: CNOT (CX) or ECR

All other gates (Hadamard, Toffoli, etc.) must be decomposed into this native set. The decomposition choice affects both circuit depth and fidelity.

# Example: Decomposing a Toffoli (CCX) gate
from qiskit import QuantumCircuit

qc = QuantumCircuit(3)
qc.ccx(0, 1, 2)  # Toffoli gate
qc_decomposed = transpile(qc, basis_gates=['sx', 'rz', 'cx'])
print(f"Decomposed Toffoli: {qc_decomposed.count_ops()}")
# Output: {'cx': 6, 'rz': 7, 'sx': 9} — 6 CNOT gates!

Dynamical Decoupling

Idle qubits accumulate decoherence errors. Dynamical decoupling inserts sequences of gate pairs (e.g., $X$-$X$ or $X$-$Y$-$X$-$Y$) during idle periods to suppress decoherence:

from qiskit import transpile
from qiskit.transpiler.passes import DynamicalDecoupling

# Apply dynamical decoupling with XX sequence
qc_dd = transpile(
    qc,
    backend=backend,
    scheduling_method='asap',
    instruction_durations=backend.instruction_durations,
)

The effect can be significant: on IBM hardware, dynamical decoupling can improve circuit fidelity by 10-30% for circuits with idle qubits.

33.13 Cross-Algorithm Comparison: Lessons from the Portfolio

Having implemented all eight algorithms, let's step back and draw comparative lessons:

Algorithm Complexity vs. Hardware Feasibility:

Key insight: The algorithms that are most feasible on NISQ hardware are either very shallow (Bell, DJ, BV) or variational (VQE, QAOA). Algorithms that require deep circuits (Shor's, QFT for large n) are not yet feasible.

33.14 Detailed Noise Models for Hardware Simulation

Understanding and modeling noise is essential for interpreting hardware results. This section provides detailed noise models that approximate real quantum hardware behavior.

Comprehensive Noise Model for IBM Quantum Processors:

from qiskit_aer.noise import (
    NoiseModel, depolarizing_error, thermal_relaxation_error,
    readout_error, phase_amplitude_damping_error
)

def create_ibm_noise_model(n_qubits=2, t1_us=300, t2_us=200,
                            cx_error=0.007, sx_error=0.0002,
                            x_error=0.0002, measure_error=0.015):
    """
    Create a realistic noise model approximating an IBM Quantum processor.

    Args:
        n_qubits: number of qubits
        t1_us: T1 time in microseconds
        t2_us: T2 time in microseconds (must be ≤ 2*T1)
        cx_error: CNOT gate error rate
        sx_error: SX gate error rate
        x_error: X gate error rate
        measure_error: readout error rate

    Returns:
        NoiseModel for Qiskit Aer
    """
    noise_model = NoiseModel()

    # Gate times (in microseconds)
    sx_time = 0.032  # 32 ns
    cx_time = 0.320  # 320 ns
    measure_time = 2.0  # 2 μs

    # Single-qubit gate errors (depolarizing + thermal relaxation)
    for gate_name, error_rate, gate_time in [
        ('sx', sx_error, sx_time),
        ('x', x_error, sx_time),
    ]:
        # Depolarizing error
        dep_err = depolarizing_error(error_rate, 1)
        # Thermal relaxation
        thermal_err = thermal_relaxation_error(t1_us * 1e3, t2_us * 1e3, 
                                                gate_time * 1e3)
        # Combined error
        combined_err = dep_err.compose(thermal_err)
        noise_model.add_all_qubit_quantum_error(combined_err, [gate_name])

    # Two-qubit gate errors
    cx_dep = depolarizing_error(cx_error, 2)
    # Thermal relaxation on both qubits during CX
    cx_thermal_0 = thermal_relaxation_error(t1_us * 1e3, t2_us * 1e3,
                                             cx_time * 1e3)
    cx_thermal_1 = thermal_relaxation_error(t1_us * 1e3, t2_us * 1e3,
                                             cx_time * 1e3)
    cx_thermal = cx_thermal_0.tensor(cx_thermal_1)
    cx_combined = cx_dep.compose(cx_thermal)
    noise_model.add_all_qubit_quantum_error(cx_combined, ['cx'])

    # Readout errors
    readout_err = readout_error([
        [1 - measure_error, measure_error],
        [measure_error, 1 - measure_error]
    ])
    noise_model.add_all_qubit_readout_error(readout_err)

    return noise_model

# Test the noise model
noise_model = create_ibm_noise_model()

# Compare ideal and noisy simulation for Bell state
from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator

qc = QuantumCircuit(2, 2)
qc.h(0)
qc.cx(0, 1)
qc.measure([0, 1], [0, 1])

# Ideal simulation
sim_ideal = AerSimulator()
result_ideal = sim_ideal.run(qc, shots=8192).result()
counts_ideal = result_ideal.get_counts()

# Noisy simulation
sim_noisy = AerSimulator(noise_model=noise_model)
result_noisy = sim_noisy.run(qc, shots=8192).result()
counts_noisy = result_noisy.get_counts()

print("Ideal results:", counts_ideal)
print("Noisy results:", counts_noisy)

# Compute fidelity
p_ideal_00 = counts_ideal.get('00', 0) / 8192
p_ideal_11 = counts_ideal.get('11', 0) / 8192
p_noisy_00 = counts_noisy.get('00', 0) / 8192
p_noisy_11 = counts_noisy.get('11', 0) / 8192
p_noisy_01 = counts_noisy.get('01', 0) / 8192
p_noisy_10 = counts_noisy.get('10', 0) / 8192

fidelity = p_noisy_00 + p_noisy_11
print(f"\nBell state fidelity: {fidelity:.4f}")
print(f"Leakage to |01>: {p_noisy_01:.4f}")
print(f"Leakage to |10>: {p_noisy_10:.4f}")

Understanding the noise model:

The realistic noise model includes three types of errors:

1. Depolarizing errors: With probability $p$, the gate output is completely random (the identity channel becomes the depolarizing channel). This models systematic gate imperfections.

2. Thermal relaxation: Energy relaxation (T1) and dephasing (T2) cause the qubit state to decay toward $|0\rangle$ and lose phase information, respectively. The amplitude damping channel models T1 decay:

$$\mathcal{E}_{\text{T1}}(\rho) = \begin{pmatrix} \rho_{00} + (1-e^{-t/T_1})\rho_{11} & e^{-t/T_2}\rho_{01} \\ e^{-t/T_2}\rho_{10}^* & e^{-t/T_1}\rho_{11} \end{pmatrix}$$

3. Readout errors: Measurement misclassification where $|0\rangle$ is reported as $|1\rangle$ (and vice versa) with probability $p_{\text{meas}}$.

These three error sources combine multiplicatively: the overall circuit fidelity is approximately the product of all individual fidelities.

33.15 Advanced Error Mitigation Techniques

Beyond the basic zero-noise extrapolation introduced in Project 7, several advanced error mitigation techniques are essential for extracting useful results from NISQ hardware.

Measurement Error Mitigation

Readout errors are often the dominant error source, especially for circuits with few gates. The calibration matrix $M$ characterizes the readout errors:

$$M = \begin{pmatrix} P(0|0) & P(0|1) \\ P(1|0) & P(1|1) \end{pmatrix}$$

where $P(j|i)$ is the probability of measuring $j$ when the true state is $|i\rangle$. For a well-calibrated system, $P(0|0) \approx 0.985$ and $P(1|1) \approx 0.980$.

The corrected counts are obtained by inverting the calibration matrix:

$$\vec{p}_{\text{corrected}} = M^{-1} \vec{p}_{\text{measured}}$$

from qiskit.result import LocalReadoutMitigator
from qiskit_aer.noise import NoiseModel, ReadoutError

# Create a readout error model
p0g1 = 0.02  # P(measure 0 | prepared 1) = 2%
p1g0 = 0.015  # P(measure 1 | prepared 0) = 1.5%

readout_err = ReadoutError([
    [1 - p1g0, p1g0],
    [p0g1, 1 - p0g1]
])

noise_model = NoiseModel()
noise_model.add_all_qubit_readout_error(readout_err)

# Run with and without mitigation
sim = AerSimulator(noise_model=noise_model)
result = sim.run(qc_bell, shots=8192).result()
counts_raw = result.get_counts()

# Apply readout mitigation manually
# For a 2-qubit system, the calibration matrix is 4x4
M = np.array([
    [(1-p1g0)**2, (1-p1g0)*p0g1, p0g1*(1-p1g0), p0g1**2],
    [(1-p1g0)*p1g0, (1-p1g0)*(1-p0g1), p0g1*p1g0, p0g1*(1-p0g1)],
    [p1g0*(1-p1g0), p1g0*p0g1, (1-p0g1)*(1-p1g0), (1-p0g1)*p0g1],
    [p1g0**2, p1g0*(1-p0g1), (1-p0g1)*p1g0, (1-p0g1)**2]
])

p_raw = np.array([counts_raw.get(f'{i:02b}', 0) for i in range(4)]) / 8192
p_corrected = np.linalg.solve(M, p_raw)

print("Raw probabilities:", dict(zip(['00','01','10','11'], p_raw.round(4))))
print("Corrected probabilities:", dict(zip(['00','01','10','11'], np.clip(p_corrected, 0, 1).round(4))))
print(f"Ideal probabilities: {{'00': 0.5, '01': 0.0, '10': 0.0, '11': 0.5}}")

Probabilistic Error Cancellation (PEC)

PEC is a more powerful but more expensive technique that constructs an unbiased estimator of the ideal expectation value by sampling from a set of "negative" and "positive" noise operations:

$$\langle O \rangle_{\text{ideal}} = \sum_i \alpha_i \langle O \rangle_{\mathcal{N}_i}$$

where $\alpha_i$ are signed weights and $\mathcal{N}_i$ are noise operations. The sampling overhead scales as $\gamma^d$ where $\gamma > 1$ is the "one-norm factor" and $d$ is the circuit depth. For typical noise levels, $\gamma \approx 5-10$, meaning PEC becomes impractical for deep circuits.

Twirled Readout Error Extinction (TREX)

A simpler approach that applies random bit-flips to the measurement basis and averages over them to cancel readout bias. This is effective for circuits where readout error is the dominant source of noise.

When to use each technique:

The practical recommendation: start with readout mitigation (cheap and effective), add zero-noise extrapolation for medium-depth circuits, and reserve PEC for short circuits where high accuracy is required.

The Variational Advantage:

VQE and QAOA are uniquely suited to NISQ hardware because:

1. Shallow circuits: Each iteration uses a shallow circuit, reducing noise impact.
2. Classical optimization: The classical optimizer can partially compensate for systematic noise.
3. Variational bound: For VQE, the estimated energy is always an upper bound on the true energy (for a reasonable ansatz), providing a built-in error check.
4. Adaptability: The ansatz can be tailored to the specific problem and hardware topology.

Example 33.5: Noise Resilience of Variational Algorithms

Consider VQE on a noisy device. The energy landscape is shifted by noise:

$$E_{\text{noisy}}(\theta) = E_{\text{ideal}}(\theta) + \epsilon_{\text{noise}}(\theta)$$

For depolarizing noise with rate $p$, the noisy energy is:

$$E_{\text{noisy}} = (1 - p)^d \cdot E_{\text{ideal}} + (1 - (1-p)^d) \cdot \text{Tr}(H) / 2^n$$

where $d$ is the circuit depth and $\text{Tr}(H) / 2^n$ is the average energy (for depolarizing noise). The key property: the noisy energy is a convex combination of the ideal energy and the average energy, so:

$$E_{\text{noisy}} \geq E_{\text{ideal}} \cdot (1 - p)^d + E_{\text{avg}} \cdot (1 - (1-p)^d)$$

If $E_{\text{avg}} > E_0$ (the ground state energy), which is typically the case, then the noisy estimate is still above the true ground state energy. This means VQE's variational principle provides a floor on the error — it can't "overshoot" and find an energy below the true ground state.

This is why VQE works on NISQ hardware: the variational principle provides error resilience that other algorithms lack.

33.16 Building Your Quantum Portfolio

The algorithms in this chapter form the foundation of a quantum computing portfolio — a collection of implemented and tested quantum programs that demonstrates your practical skills. Here's how to organize and extend this portfolio:

Portfolio Structure:

quantum-portfolio/
├── 01-bell-state/
│   ├── bell_state.ipynb          # Implementation and analysis
│   ├── bell_state_results.json    # Hardware results
│   └── README.md                  # Summary and key findings
├── 02-deutsch-jozsa/
│   ├── deutsch_jozsa.ipynb
│   ├── dj_results.json
│   └── README.md
├── 03-bernstein-vazirani/
│   ├── bv.ipynb
│   ├── bv_results.json
│   └── README.md
├── 04-qft/
│   ├── qft.ipynb
│   ├── qft_results.json
│   └── README.md
├── 05-grover/
│   ├── grover.ipynb
│   ├── grover_results.json
│   └── README.md
├── 06-shor/
│   ├── shor_15.ipynb
│   ├── shor_results.json
│   └── README.md
├── 07-vqe/
│   ├── vqe_h2.ipynb
│   ├── vqe_results.json
│   └── README.md
├── 08-qaoa/
│   ├── qaoa_maxcut.ipynb
│   ├── qaoa_results.json
│   └── README.md
└── README.md                      # Portfolio overview

For each algorithm, include: - Mathematical derivation (in LaTeX) - Circuit diagram (ASCII art) - Complete Qiskit implementation - Ideal simulation results - Hardware results (with backend name and date) - Noise analysis and fidelity comparison - Discussion of discrepancies between ideal and hardware results

Extending the portfolio: - Add new algorithms as you learn them (quantum phase estimation, quantum walk, quantum key distribution) - Compare results across different hardware backends - Implement error mitigation techniques and measure their impact - Contribute your best implementations to open-source projects