Qiskit Beginner Free 15/61 in series 35 minutes

Debugging Quantum Circuits with Qiskit's Statevector Simulator

Use the statevector simulator to debug quantum circuits: print amplitudes, plot probability distributions, probe mid-circuit state, and catch common bugs.

What you'll learn

  • Qiskit
  • debugging
  • statevector
  • circuit visualization
  • probabilities
  • intermediate measurement

Prerequisites

  • Basic Python (variables, functions, loops)
  • No quantum physics background needed

Why Quantum Debugging Is Hard

Debugging quantum circuits is fundamentally different from debugging classical code. Before diving into techniques, it helps to understand why quantum bugs are so tricky to find and fix.

No print-statement debugging

In classical programming, you can insert print() statements anywhere to inspect variables. Quantum mechanics forbids this. Measuring a quantum state collapses it, destroying the superposition you want to examine. If you add a measurement to debug a circuit, you change the circuit’s behavior.

This is why simulators are essential for debugging. The statevector simulator sidesteps the measurement problem entirely by computing the full quantum state mathematically, without any physical collapse.

Exponential state space

A classical register of n bits has exactly n values to inspect. A quantum register of n qubits has 2^n complex amplitudes. For 10 qubits, that is 1,024 amplitudes. For 20 qubits, over one million. For 30 qubits, over one billion. You cannot scroll through a billion complex numbers and spot the bug by eye. Effective quantum debugging requires targeted checks, not manual inspection.

Phase errors are invisible to measurements

Two quantum states can produce identical measurement probability distributions yet behave completely differently when used as part of a larger algorithm. This happens when the states differ only in their phases. Consider this concrete example:

from qiskit import QuantumCircuit
from qiskit.quantum_info import Statevector, state_fidelity
import numpy as np

# Correct circuit: prepare |+> state
qc_correct = QuantumCircuit(1)
qc_correct.h(0)

# Bugged circuit: prepare |-> state (has a phase error)
qc_bugged = QuantumCircuit(1)
qc_bugged.x(0)
qc_bugged.h(0)

sv_correct = Statevector(qc_correct)
sv_bugged = Statevector(qc_bugged)

print("Correct |+> amplitudes:", sv_correct.data)
# Output: [0.70710678+0.j, 0.70710678+0.j]
print("Bugged  |-> amplitudes:", sv_bugged.data)
# Output: [0.70710678+0.j, -0.70710678+0.j]

print("\nProbabilities match:",
      np.allclose(sv_correct.probabilities(), sv_bugged.probabilities()))
# Output: True

# Now use each as input to another Hadamard (interference step)
qc_correct.h(0)
qc_bugged.h(0)

sv_after_correct = Statevector(qc_correct)
sv_after_bugged = Statevector(qc_bugged)

print("\nAfter interference (H gate):")
print("Correct final state:", sv_after_correct.data)
# Output: [1.+0.j, 0.+0.j]  -> measures |0> with 100% probability
print("Bugged final state: ", sv_after_bugged.data)
# Output: [0.+0.j, 1.+0.j]  -> measures |1> with 100% probability

The two circuits produce identical histograms before the interference step. After the interference step, they give completely opposite results. This is exactly what happens in algorithms like Grover’s search or the Quantum Fourier Transform, where phase errors silently propagate and cause wrong answers at the end.

Non-determinism in measurement

A single measurement of a quantum state gives one bitstring sampled from the probability distribution. Running 100 shots gives a noisy approximation. You need thousands of shots to distinguish a correct distribution from a slightly bugged one. Statevector debugging avoids this entirely by giving you the exact state.

Setup

pip install qiskit qiskit-aer matplotlib
from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator
from qiskit.quantum_info import Statevector, state_fidelity
from qiskit.visualization import plot_bloch_multivector, plot_state_city
import numpy as np

Technique 1: Print State Amplitudes

The simplest debugging tool is printing the full statevector after a circuit runs:

# Build a Bell state
qc = QuantumCircuit(2)
qc.h(0)
qc.cx(0, 1)

sv = Statevector(qc)
print("Statevector:")
print(sv)
print()
print("Amplitudes:")
for i, amp in enumerate(sv):
    bits = format(i, f"0{qc.num_qubits}b")
    print(f"  |{bits}> : {amp:.4f}")

For a correct Bell state you expect:

  • |00> amplitude: 0.7071 + 0j
  • |11> amplitude: 0.7071 + 0j
  • |01> and |10>: 0

If any of these differ, the circuit has a bug.

Technique 2: Probability Distributions

Amplitudes include phase information, but measurements only see probabilities. To check if your circuit will produce the right measurement statistics:

probs = sv.probabilities_dict()
print("Measurement probabilities:")
for bitstring, prob in sorted(probs.items()):
    bar = "#" * int(prob * 40)
    print(f"  |{bitstring}> : {prob:.4f}  {bar}")

Plot the probability distribution visually:

from qiskit.visualization import plot_distribution

# Simulate with shots to compare to ideal
sim = AerSimulator()
qc_m = qc.copy()
qc_m.measure_all()
job = sim.run(qc_m, shots=2000)
counts = job.result().get_counts()

fig = plot_distribution(counts)
fig.suptitle("Bell State Measurement Distribution")
fig.savefig("bell_distribution.png")
print("Saved bell_distribution.png")

Technique 3: Mid-Circuit Probes with save_statevector()

The most powerful debugging technique inserts save_statevector() instructions at specific points in the circuit. This is equivalent to putting a breakpoint inside a classical program:

qc = QuantumCircuit(2)

qc.h(0)
qc.save_statevector(label="after_h")  # probe 1

qc.cx(0, 1)
qc.save_statevector(label="after_cx")  # probe 2

qc.x(1)
qc.save_statevector(label="after_x")  # probe 3

sim = AerSimulator()
result = sim.run(qc).result()

for label in ["after_h", "after_cx", "after_x"]:
    sv = result.data(0)[label]
    print(f"\nState {label}:")
    for i, amp in enumerate(sv):
        bits = format(i, f"0{qc.num_qubits}b")
        if abs(amp) > 1e-10:
            print(f"  |{bits}> : {amp:.4f}")

This reveals exactly how the state evolves step by step. If after_cx looks wrong, the bug is in the H gate or the initial state, not in subsequent gates.

Technique 4: Checking for Expected Bell and GHZ States

For common entangled states, verify against the known exact statevector:

from qiskit.quantum_info import Statevector, state_fidelity

def check_bell_state(qc: QuantumCircuit, which: str = "phi_plus") -> float:
    """Compute fidelity against a target Bell state."""
    bell_states = {
        "phi_plus":  Statevector([1/np.sqrt(2), 0, 0,  1/np.sqrt(2)]),
        "phi_minus": Statevector([1/np.sqrt(2), 0, 0, -1/np.sqrt(2)]),
        "psi_plus":  Statevector([0, 1/np.sqrt(2),  1/np.sqrt(2), 0]),
        "psi_minus": Statevector([0, 1/np.sqrt(2), -1/np.sqrt(2), 0]),
    }
    target = bell_states[which]
    actual = Statevector(qc)
    return state_fidelity(actual, target)

# Correct Bell circuit
qc_correct = QuantumCircuit(2)
qc_correct.h(0)
qc_correct.cx(0, 1)

# Bugged Bell circuit (wrong qubit order in CX)
qc_bugged = QuantumCircuit(2)
qc_bugged.h(0)
qc_bugged.cx(1, 0)  # control and target swapped

print(f"Correct circuit fidelity: {check_bell_state(qc_correct):.6f}")
# Output: Correct circuit fidelity: 1.000000
print(f"Bugged circuit fidelity:  {check_bell_state(qc_bugged):.6f}")
# Output: Bugged circuit fidelity:  0.250000

Fidelity of 1.0 means the circuit is correct. Anything less indicates a bug.

Technique 5: Bloch Sphere Visualization

For single-qubit states, the Bloch sphere gives an immediate geometric intuition:

# Trace over one qubit to get single-qubit Bloch sphere
qc = QuantumCircuit(1)
qc.h(0)
qc.s(0)   # S gate: phase rotation by pi/2

sv = Statevector(qc)
fig = plot_bloch_multivector(sv)
fig.savefig("bloch.png")
print("State on Bloch sphere saved to bloch.png")
print(f"Amplitudes: {sv.data}")
# Output: Amplitudes: [0.70710678+0.j 0.        +0.70710678j]

Bloch sphere positions to know:

  • |0> is the north pole: (0, 0, 1)
  • |+> (after H) is on the equator: (1, 0, 0)
  • |+i> (after H then S) is on the equator: (0, 1, 0)
  • Any off-axis position reveals an unexpected phase

Technique 6: City and Hinton Plots

For multi-qubit states, city and Hinton plots show the real and imaginary parts of the density matrix:

from qiskit.quantum_info import DensityMatrix
from qiskit.visualization import plot_state_city, plot_state_hinton

# GHZ state
qc_ghz = QuantumCircuit(3)
qc_ghz.h(0)
qc_ghz.cx(0, 1)
qc_ghz.cx(1, 2)

dm = DensityMatrix(qc_ghz)
fig_city   = plot_state_city(dm, title="GHZ State - City Plot")
fig_hinton = plot_state_hinton(dm, title="GHZ State - Hinton Plot")
fig_city.savefig("ghz_city.png")
fig_hinton.savefig("ghz_hinton.png")

The city plot shows amplitude magnitudes as bar heights. For a GHZ state you expect four tall bars at the corners corresponding to |000><000|, |000><111|, |111><000|, and |111><111|. If bars appear elsewhere, the circuit has an error.

Statevector vs Density Matrix

When should you use Statevector and when should you use DensityMatrix? The answer depends on whether your circuit involves noise or mid-circuit measurements.

Statevector represents a pure quantum state. It works for any circuit that consists entirely of unitary gates (no noise, no classical conditioning, no mid-circuit measurement). It stores 2^n complex amplitudes.

DensityMatrix represents both pure and mixed quantum states. A mixed state arises when there is classical uncertainty about which quantum state the system is in, which happens after noise or after tracing out part of a system. It stores a 2^n by 2^n complex matrix.

Comparing Statevector and DensityMatrix for a Bell state

from qiskit import QuantumCircuit
from qiskit.quantum_info import Statevector, DensityMatrix
import numpy as np

# Bell state circuit
qc = QuantumCircuit(2)
qc.h(0)
qc.cx(0, 1)

# Statevector representation
sv = Statevector(qc)
print("Statevector:", sv.data)
# Output: [0.70710678+0.j, 0.+0.j, 0.+0.j, 0.70710678+0.j]

# DensityMatrix representation
dm = DensityMatrix(qc)
print("\nDensity matrix:\n", np.round(dm.data, 4))
# Output: 4x4 matrix with 0.5 at corners (|00><00|, |00><11|, |11><00|, |11><11|)

For a pure state, the density matrix equals the outer product of the statevector with itself: rho = |psi><psi|. Both representations contain the same information.

When noise makes Statevector misleading

When you add noise to a circuit, the state becomes mixed. Statevector cannot represent mixed states and gives incorrect results. DensityMatrix handles noise correctly.

from qiskit.quantum_info import DensityMatrix, Statevector
from qiskit_aer.noise import depolarizing_error
from qiskit import QuantumCircuit
import numpy as np

# Start with a simple |+> state
qc = QuantumCircuit(1)
qc.h(0)

# Pure state: Statevector and DensityMatrix agree
sv_pure = Statevector(qc)
dm_pure = DensityMatrix(qc)
print("Pure state purity:", dm_pure.purity())
# Output: 1.0

# Now simulate the effect of depolarizing noise on the density matrix.
# A depolarizing channel with parameter p maps rho to (1-p)*rho + p*I/2.
p_noise = 0.3
dm_noisy = DensityMatrix(
    (1 - p_noise) * dm_pure.data + p_noise * np.eye(2) / 2
)

print("Noisy state purity:", dm_noisy.purity())
# Output: approximately 0.755 (less than 1.0, so the state is mixed)

print("\nNoisy density matrix diagonal:", np.diag(dm_noisy.data).real)
# The diagonal entries show probabilities are no longer [0.5, 0.5]
# but are slightly shifted toward the maximally mixed state

Use the purity() method to check whether a state is pure or mixed. A purity of 1.0 means a pure state. Anything less than 1.0 indicates a mixed state, and you should use DensityMatrix for accurate debugging.

Debugging Qubit Ordering Confusion

Qiskit uses little-endian qubit ordering. This is the single most common source of confusion and bugs for new Qiskit users. Understanding it thoroughly saves hours of debugging.

The rule is: qubit 0 is the rightmost bit in a state label string. A 2-qubit state where qubit 0 is |1> and qubit 1 is |0> is written as '01' in Qiskit output, because qubit 1 (leftmost position) is 0 and qubit 0 (rightmost position) is 1.

Demonstrating qubit ordering

from qiskit import QuantumCircuit
from qiskit.quantum_info import Statevector

# Apply X only to qubit 0
qc = QuantumCircuit(2)
qc.x(0)

sv = Statevector(qc)
print("After X on qubit 0:")
for i, amp in enumerate(sv):
    bits = format(i, f"0{qc.num_qubits}b")
    if abs(amp) > 1e-10:
        print(f"  |{bits}> has amplitude {amp:.4f}")
# Output: |01> has amplitude 1.0000
# The '1' is in the rightmost position, confirming qubit 0 is rightmost

# Apply X only to qubit 1
qc2 = QuantumCircuit(2)
qc2.x(1)

sv2 = Statevector(qc2)
print("\nAfter X on qubit 1:")
for i, amp in enumerate(sv2):
    bits = format(i, f"0{qc2.num_qubits}b")
    if abs(amp) > 1e-10:
        print(f"  |{bits}> has amplitude {amp:.4f}")
# Output: |10> has amplitude 1.0000
# The '1' is in the leftmost position, confirming qubit 1 is leftmost

Interpreting 3-qubit measurement strings

For a 3-qubit system, the bitstring '101' means:

  • Qubit 2 (leftmost) = 1
  • Qubit 1 (middle) = 0
  • Qubit 0 (rightmost) = 1
qc3 = QuantumCircuit(3)
qc3.x(0)  # Set qubit 0 to |1>
qc3.x(2)  # Set qubit 2 to |1>
# Qubit 1 stays |0>

sv3 = Statevector(qc3)
probs = sv3.probabilities_dict()
print("Measurement probabilities:", probs)
# Output: {'101': 1.0}
# Reading right to left: qubit 0=1, qubit 1=0, qubit 2=1

Helper function for labeled amplitude printing

This utility function prints amplitudes with explicit qubit labels to eliminate ordering confusion:

from qiskit.quantum_info import Statevector

def print_labeled_state(sv: Statevector, qubit_names: list = None):
    """Print statevector amplitudes with explicit qubit labels.

    Args:
        sv: The statevector to print.
        qubit_names: Optional names for each qubit (index 0 = qubit 0).
                     Defaults to q0, q1, q2, etc.
    """
    n = int(np.log2(len(sv)))
    if qubit_names is None:
        qubit_names = [f"q{i}" for i in range(n)]

    print(f"{'Basis':>10}  {'Amplitude':>20}  {'Prob':>8}  Qubit values")
    print("-" * 70)
    for i, amp in enumerate(sv):
        if abs(amp) < 1e-10:
            continue
        bits = format(i, f"0{n}b")
        prob = abs(amp) ** 2
        # Map each bit to its qubit label (remember: reversed order)
        labels = [f"{qubit_names[j]}={bits[n - 1 - j]}"
                  for j in range(n)]
        print(f"  |{bits}>  {amp:>20.4f}  {prob:>8.4f}  {', '.join(labels)}")

# Example: Bell state
qc_bell = QuantumCircuit(2)
qc_bell.h(0)
qc_bell.cx(0, 1)
sv_bell = Statevector(qc_bell)
print_labeled_state(sv_bell, qubit_names=["control", "target"])

Systematic Gate Debugging Workflow

When a circuit produces wrong results, resist the urge to stare at the code and guess. Follow this systematic 5-step workflow instead.

The 5 steps

  1. Verify the input state. Check that all qubits start in the expected state.
  2. Test each gate in isolation. Apply each gate to a known input and verify the output.
  3. Add statevector probes at key points. Insert save_statevector after each logical section of the circuit.
  4. Compare against the expected analytical result. Compute the expected statevector by hand or with a reference implementation.
  5. Run a fidelity check. Use state_fidelity to get a single pass/fail number.

Applying the workflow to a buggy GHZ circuit

A 3-qubit GHZ state is (|000> + |111>) / sqrt(2). The correct circuit applies H to qubit 0, then CNOT from 0 to 1, then CNOT from 1 to 2. Here is a buggy version with the CNOTs in the wrong order:

from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator
from qiskit.quantum_info import Statevector, state_fidelity
import numpy as np

# Buggy GHZ circuit: CNOT order is wrong
qc_bug = QuantumCircuit(3)
qc_bug.h(0)
qc_bug.cx(0, 2)  # Bug: should be cx(0, 1) first
qc_bug.cx(1, 2)  # Bug: should be cx(1, 2) but qubit 1 is not entangled yet

# Step 1: Verify input state
sv_init = Statevector.from_label("000")
print("Step 1 - Input state:", sv_init.data)
# Output: [1, 0, 0, 0, 0, 0, 0, 0]  (correct: all qubits in |0>)

# Step 2: Test H gate in isolation
qc_h_only = QuantumCircuit(3)
qc_h_only.h(0)
sv_after_h = Statevector(qc_h_only)
expected_after_h = Statevector([1/np.sqrt(2), 1/np.sqrt(2), 0, 0, 0, 0, 0, 0])
print(f"\nStep 2 - H gate fidelity: {state_fidelity(sv_after_h, expected_after_h):.6f}")
# Output: 1.000000 (H gate works correctly)

# Step 3: Add probes at key points
qc_probed = QuantumCircuit(3)
qc_probed.h(0)
qc_probed.save_statevector(label="after_h")
qc_probed.cx(0, 2)
qc_probed.save_statevector(label="after_cx1")
qc_probed.cx(1, 2)
qc_probed.save_statevector(label="after_cx2")

sim = AerSimulator()
result = sim.run(qc_probed).result()

print("\nStep 3 - Probing intermediate states:")
for label in ["after_h", "after_cx1", "after_cx2"]:
    sv_probe = result.data(0)[label]
    print(f"\n  {label}:")
    for i, amp in enumerate(sv_probe):
        if abs(amp) > 1e-10:
            bits = format(i, "03b")
            print(f"    |{bits}> : {amp:.4f}")

# Step 4: Compare against expected GHZ state
# After after_cx1, we expect (|000> + |101>)/sqrt(2) instead of (|000> + |110>)/sqrt(2)
# This tells us the first CNOT targets qubit 2 instead of qubit 1

# Step 5: Fidelity check
target_ghz = Statevector([1/np.sqrt(2), 0, 0, 0, 0, 0, 0, 1/np.sqrt(2)])
sv_final = Statevector(qc_bug)
fidelity = state_fidelity(sv_final, target_ghz)
print(f"\nStep 5 - Final fidelity against GHZ: {fidelity:.6f}")
# Output: less than 1.0, confirming the bug

The probes at step 3 reveal the bug clearly. After the first CNOT, the state is (|000> + |101>)/sqrt(2) instead of the expected (|000> + |110>)/sqrt(2). This tells you that the first CNOT targets qubit 2 when it should target qubit 1. Fixing the CNOT order to cx(0,1) then cx(1,2) produces the correct GHZ state.

The save_statevector technique becomes especially valuable in multi-step algorithms where you need to verify that each stage works correctly. Here we implement Grover’s algorithm for 3 qubits with a single marked state and probe the state after each key step.

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

def grover_oracle(qc, marked_state, n):
    """Apply an oracle that flips the phase of the marked state."""
    # Convert marked state to binary and apply X gates to zero bits
    marked_bits = format(marked_state, f"0{n}b")
    for i in range(n):
        if marked_bits[n - 1 - i] == '0':
            qc.x(i)
    # Multi-controlled Z gate: apply H to last qubit, MCX, then H again
    qc.h(n - 1)
    qc.mcx(list(range(n - 1)), n - 1)
    qc.h(n - 1)
    # Undo X gates
    for i in range(n):
        if marked_bits[n - 1 - i] == '0':
            qc.x(i)

def grover_diffusion(qc, n):
    """Apply the Grover diffusion operator."""
    for i in range(n):
        qc.h(i)
    for i in range(n):
        qc.x(i)
    qc.h(n - 1)
    qc.mcx(list(range(n - 1)), n - 1)
    qc.h(n - 1)
    for i in range(n):
        qc.x(i)
    for i in range(n):
        qc.h(i)

n = 3
marked_state = 5  # Binary: 101

# Build Grover circuit with probes
qc = QuantumCircuit(n)

# Initialize: apply H to all qubits
for i in range(n):
    qc.h(i)
qc.save_statevector(label="after_init")

# Grover iteration 1
grover_oracle(qc, marked_state, n)
qc.save_statevector(label="after_oracle_1")

grover_diffusion(qc, n)
qc.save_statevector(label="after_diffusion_1")

# Grover iteration 2
grover_oracle(qc, marked_state, n)
qc.save_statevector(label="after_oracle_2")

grover_diffusion(qc, n)
qc.save_statevector(label="after_diffusion_2")

# Run simulation
sim = AerSimulator()
result = sim.run(qc).result()

# Check marked state amplitude at each probe point
marked_bits = format(marked_state, f"0{n}b")
print(f"Tracking amplitude of marked state |{marked_bits}> (state index {marked_state})")
print()

for label in ["after_init", "after_oracle_1", "after_diffusion_1",
              "after_oracle_2", "after_diffusion_2"]:
    sv = result.data(0)[label]
    amp = sv[marked_state]
    prob = abs(amp) ** 2
    print(f"  {label:>25}: amplitude = {amp:.4f}, probability = {prob:.4f}")

Expected output shows the marked state’s probability growing after iteration 1 and being close to maximum after iteration 1 (for n=3, the optimal number of iterations is approximately pi/4 * sqrt(8) ~ 2.2, so 2 iterations is near optimal). After too many iterations, the amplitude starts shrinking. This probing approach lets you diagnose whether the oracle or diffusion operator is functioning correctly and whether you are using the right number of iterations.

For n=3 with 1 marked state:

  • After initialization: probability = 1/8 = 0.125 (uniform superposition)
  • After oracle 1: probability = 0.125 (oracle only flips phase, does not change probability)
  • After diffusion 1: probability ~ 0.78 (amplitude amplification working)
  • After oracle 2: probability ~ 0.78 (again, oracle preserves probability)
  • After diffusion 2: probability ~ 0.95 (near maximum)

If you see that the probability does not increase after the diffusion step, the diffusion operator has a bug. If the probability does not stay constant after the oracle step, the oracle is changing amplitudes instead of just phases.

Tensor Product Structure and Factorizability

When debugging quantum circuits, it helps to verify whether a multi-qubit state is entangled or factorizable. A product state like |0>|+> factors into individual qubit states. An entangled state like the Bell state does not factor. If your circuit should produce entanglement but the state factors, or vice versa, you have a bug.

Detecting entanglement with partial trace and entropy

The von Neumann entropy of a subsystem tells you whether it is entangled with the rest of the system. For a product state, the reduced density matrix of each subsystem is pure (entropy = 0). For an entangled state, the reduced density matrix is mixed (entropy > 0).

from qiskit import QuantumCircuit
from qiskit.quantum_info import Statevector, DensityMatrix, partial_trace, entropy
import numpy as np

def check_entanglement(qc, qubit_a=0, qubit_b=1):
    """Check if two qubits are entangled by computing subsystem entropy."""
    dm = DensityMatrix(qc)
    n = qc.num_qubits
    # Trace out all qubits except qubit_a
    qubits_to_trace = [i for i in range(n) if i != qubit_a]
    reduced_dm = partial_trace(dm, qubits_to_trace)
    ent = entropy(reduced_dm, base=2)
    return ent

# Test 1: Product state |00> (not entangled)
qc_product = QuantumCircuit(2)
# No gates: state is |00>
ent_product = check_entanglement(qc_product)
print(f"|00> state entropy: {ent_product:.4f}")
# Output: 0.0000 (product state, no entanglement)

# Test 2: Bell state (maximally entangled)
qc_bell = QuantumCircuit(2)
qc_bell.h(0)
qc_bell.cx(0, 1)
ent_bell = check_entanglement(qc_bell)
print(f"Bell state entropy: {ent_bell:.4f}")
# Output: 1.0000 (maximally entangled for 2 qubits)

# Test 3: |+>|+> state (looks complex but is NOT entangled)
qc_plus_plus = QuantumCircuit(2)
qc_plus_plus.h(0)
qc_plus_plus.h(1)
ent_plus = check_entanglement(qc_plus_plus)
print(f"|+>|+> state entropy: {ent_plus:.4f}")
# Output: 0.0000 (this is a product state despite having 4 nonzero amplitudes)

# Verify: |+>|+> = (|00> + |01> + |10> + |11>)/2
sv_pp = Statevector(qc_plus_plus)
print(f"\n|+>|+> amplitudes: {np.round(sv_pp.data, 4)}")
# All four amplitudes are 0.5. This looks like it could be entangled,
# but the entropy confirms it factors as |+> tensor |+>.

This technique is especially useful for circuits with many qubits. If your 10-qubit circuit should produce a product state between qubits 0-4 and qubits 5-9, you can verify by checking the entropy of the partition. Any nonzero entropy reveals unexpected entanglement, which points to a stray CNOT or other two-qubit gate acting across the partition.

Global Phase vs Relative Phase

When comparing statevectors, it is important to distinguish global phase from relative phase. Confusing the two leads to false bug reports or missed bugs.

Global phase is physically unobservable

A global phase is a single complex factor e^(i*theta) multiplying the entire statevector. It does not affect any measurement outcome, interference pattern, or physical observable. Two states that differ only by a global phase are physically identical.

Relative phase affects interference

A relative phase is a difference in phase between basis state amplitudes. This is physically meaningful and affects how the state interferes in subsequent operations.

from qiskit.quantum_info import Statevector, state_fidelity
import numpy as np

# State 1: |+> = (|0> + |1>)/sqrt(2)
state_plus = Statevector([1/np.sqrt(2), 1/np.sqrt(2)])

# State 2: global phase of pi/4 applied to |+>
global_phase = np.exp(1j * np.pi / 4)
state_global = Statevector(global_phase * state_plus.data)

# State 3: relative phase (|-> = (|0> - |1>)/sqrt(2))
state_minus = Statevector([1/np.sqrt(2), -1/np.sqrt(2)])

# Compare measurement probabilities
print("Probabilities of |+>:          ", state_plus.probabilities_dict())
print("Probabilities of e^(ipi/4)|+>: ", state_global.probabilities_dict())
print("Probabilities of |->:          ", state_minus.probabilities_dict())
# All three have {0: 0.5, 1: 0.5}. Probabilities cannot distinguish them.

# Compare with state_fidelity
fid_global = state_fidelity(state_plus, state_global)
fid_relative = state_fidelity(state_plus, state_minus)
print(f"\nFidelity |+> vs e^(ipi/4)|+>: {fid_global:.6f}")
# Output: 1.000000 (state_fidelity ignores global phase)
print(f"Fidelity |+> vs |->:          {fid_relative:.6f}")
# Output: 0.000000 (relative phase makes these orthogonal states)

Key takeaway: state_fidelity correctly ignores global phase differences. A fidelity of 1.0 means the states are physically identical, even if the raw amplitudes look different due to a global phase factor. A fidelity less than 1.0 means there is a real physical difference (a relative phase or amplitude difference).

Detecting relative phase bugs in practice

from qiskit import QuantumCircuit
from qiskit.quantum_info import Statevector, state_fidelity

# Correct: prepare |+i> = (|0> + i|1>)/sqrt(2) using H then S
qc_correct = QuantumCircuit(1)
qc_correct.h(0)
qc_correct.s(0)

# Bugged: accidentally used Sdg (S-dagger) instead of S
qc_bugged = QuantumCircuit(1)
qc_bugged.h(0)
qc_bugged.sdg(0)

sv_correct = Statevector(qc_correct)
sv_bugged = Statevector(qc_bugged)

print("Correct |+i>:", sv_correct.data)
# Output: [0.707+0.j, 0.+0.707j]
print("Bugged  |-i>:", sv_bugged.data)
# Output: [0.707+0.j, 0.-0.707j]

print(f"\nProbabilities match: {np.allclose(sv_correct.probabilities(), sv_bugged.probabilities())}")
# Output: True (identical measurement statistics)
print(f"Fidelity: {state_fidelity(sv_correct, sv_bugged):.6f}")
# Output: 0.000000 (states are orthogonal due to relative phase!)

These two states are orthogonal despite having identical measurement probabilities. The state_fidelity check catches the bug immediately.

Debugging Transpiled Circuits

When you run a circuit on real hardware or a noise-aware simulator, Qiskit transpiles the circuit to match the hardware’s native gate set, connectivity, and optimization preferences. This transformation can reorder qubits, decompose gates into sequences of native gates, and add SWAP gates. Verifying that transpilation preserves the intended unitary is a critical debugging step.

from qiskit import QuantumCircuit, transpile
from qiskit.quantum_info import Statevector, state_fidelity
from qiskit_aer import AerSimulator

# Original circuit: prepare a Bell state with extra single-qubit gates
qc_original = QuantumCircuit(2)
qc_original.h(0)
qc_original.cx(0, 1)
qc_original.rz(np.pi / 4, 0)
qc_original.ry(np.pi / 3, 1)

sv_original = Statevector(qc_original)
print("Original circuit statevector:")
print(sv_original)

# Transpile with different optimization levels
sim = AerSimulator()

for opt_level in [0, 1, 2, 3]:
    qc_transpiled = transpile(qc_original, backend=sim,
                               optimization_level=opt_level)
    sv_transpiled = Statevector(qc_transpiled)
    fid = state_fidelity(sv_original, sv_transpiled)
    n_gates = qc_transpiled.size()
    print(f"  optimization_level={opt_level}: "
          f"fidelity={fid:.6f}, gates={n_gates}")

All optimization levels should produce fidelity 1.0 (or extremely close to 1.0 due to floating-point precision). If fidelity drops below 0.9999, something went wrong in the transpilation process. This can happen if you pass incorrect basis gates, a faulty coupling map, or if the transpiler encounters an edge case.

Note that optimization_level=0 performs minimal transformations, making it easier to trace the mapping between original and transpiled gates. Level 3 performs aggressive optimizations that may produce a very different-looking circuit, but the final state should be identical up to global phase.

DensityMatrix Debugging with partial_trace

The partial_trace function lets you inspect the quantum state of any subset of qubits in a multi-qubit system. This is invaluable for debugging circuits where you need to verify the state of individual qubits or small groups of qubits.

from qiskit import QuantumCircuit
from qiskit.quantum_info import DensityMatrix, partial_trace, entropy, state_fidelity, Statevector
import numpy as np

# Build a 4-qubit circuit
qc = QuantumCircuit(4)
qc.h(0)           # Put qubit 0 in |+>
qc.cx(0, 1)       # Entangle qubits 0 and 1 (Bell pair)
qc.h(2)           # Put qubit 2 in |+>
qc.cx(2, 3)       # Entangle qubits 2 and 3 (Bell pair)

dm = DensityMatrix(qc)

# (a) Check if qubit 0 is in the |+> state
# Trace out qubits 1, 2, 3 to get the reduced state of qubit 0
dm_q0 = partial_trace(dm, [1, 2, 3])
target_plus = DensityMatrix.from_label("+")
fid_q0 = state_fidelity(dm_q0, target_plus)
print(f"Qubit 0 fidelity with |+>: {fid_q0:.4f}")
# Output: 0.5000
# This is NOT 1.0 because qubit 0 is entangled with qubit 1.
# Its reduced state is the maximally mixed state, not |+>.

# (b) Check if qubit 2 is unentangled from qubits 0 and 1
# Compute entropy of the subsystem {0, 1} vs {2, 3}
dm_01 = partial_trace(dm, [2, 3])
dm_23 = partial_trace(dm, [0, 1])
ent_01 = entropy(dm_01, base=2)
ent_23 = entropy(dm_23, base=2)
print(f"\nEntropy of qubits 0,1: {ent_01:.4f}")
print(f"Entropy of qubits 2,3: {ent_23:.4f}")
# Both are 1.0 because each pair is internally entangled.

# But are the two pairs entangled with EACH OTHER?
# Check by computing entropy of the {0,1} vs {2,3} partition
# If the full state is a product of the two pairs, then the
# entropy of the {0,1} subsystem should equal 1.0 (from internal entanglement)
# and the mutual information between the pairs should be zero.
dm_full = dm
ent_full = entropy(dm_full, base=2)
print(f"Entropy of full state: {ent_full:.4f}")
# Output: 0.0000 (the full state is pure)

# (c) Verify that qubit 2 is in a known state after the H gate
# by checking BEFORE the CNOT on qubits 2,3
qc_partial = QuantumCircuit(4)
qc_partial.h(0)
qc_partial.cx(0, 1)
qc_partial.h(2)
# Stop here, before cx(2, 3)

dm_partial = DensityMatrix(qc_partial)
dm_q2 = partial_trace(dm_partial, [0, 1, 3])
fid_q2_plus = state_fidelity(dm_q2, target_plus)
print(f"\nQubit 2 fidelity with |+> (before CNOT): {fid_q2_plus:.4f}")
# Output: 1.0000 (qubit 2 is in |+> and unentangled from the rest)

ent_q2 = entropy(dm_q2, base=2)
print(f"Qubit 2 entropy (before CNOT): {ent_q2:.4f}")
# Output: 0.0000 (pure state, not entangled with anything)

This technique scales to any number of qubits, though the density matrix computation itself becomes expensive beyond approximately 15 qubits.

Performance and Scaling of Debugging Tools

Statevector debugging is powerful but does not scale to large circuits. Understanding the practical limits helps you choose the right debugging strategy.

Timing and memory benchmarks

from qiskit import QuantumCircuit
from qiskit.quantum_info import Statevector
import time

for n_qubits in [5, 10, 15, 20, 25]:
    qc = QuantumCircuit(n_qubits)
    for i in range(n_qubits):
        qc.h(i)
    # Apply some entangling gates
    for i in range(n_qubits - 1):
        qc.cx(i, i + 1)

    start = time.time()
    sv = Statevector(qc)
    elapsed = time.time() - start

    n_amplitudes = 2 ** n_qubits
    memory_bytes = n_amplitudes * 16  # 16 bytes per complex128
    memory_mb = memory_bytes / (1024 * 1024)

    print(f"  n={n_qubits:>2}: {n_amplitudes:>12,} amplitudes, "
          f"{memory_mb:>10.1f} MB, {elapsed:.4f} seconds")

Approximate scaling (varies by hardware):

  • n=10: 1,024 amplitudes, ~0.02 MB, less than 1 ms
  • n=15: 32,768 amplitudes, ~0.5 MB, a few ms
  • n=20: 1,048,576 amplitudes, ~16 MB, under 1 second
  • n=25: 33,554,432 amplitudes, ~512 MB, several seconds
  • n=30: 1,073,741,824 amplitudes, ~16 GB, minutes (if you have enough RAM)

Beyond 30 qubits, statevector simulation is impractical on consumer hardware.

Alternatives for larger circuits

When your circuit exceeds the limits of statevector simulation, consider these alternatives:

Sampling-based testing. Run the circuit with many shots and compare the measurement distribution against the expected distribution using statistical tests. This does not give you the full state, but it catches most bugs.

from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator
from scipy.stats import chisquare

def test_distribution(qc, expected_probs, shots=10000, significance=0.01):
    """Test if a circuit's measurement distribution matches expected probabilities."""
    qc_m = qc.copy()
    qc_m.measure_all()
    sim = AerSimulator()
    counts = sim.run(qc_m, shots=shots).result().get_counts()

    # Build observed and expected frequency arrays
    all_bitstrings = sorted(expected_probs.keys())
    observed = [counts.get(b, 0) for b in all_bitstrings]
    expected = [expected_probs[b] * shots for b in all_bitstrings]

    chi2, p_value = chisquare(observed, f_exp=expected)
    passed = p_value > significance
    print(f"Chi-squared test: chi2={chi2:.2f}, p={p_value:.4f}, "
          f"{'PASS' if passed else 'FAIL'}")
    return passed

# Test a Bell state
qc_bell = QuantumCircuit(2)
qc_bell.h(0)
qc_bell.cx(0, 1)
test_distribution(qc_bell, {"00": 0.5, "11": 0.5})

Unit testing subsystems. Break your large circuit into smaller modules of 10 qubits or fewer. Test each module in isolation using statevector debugging. Then integrate and test the combined circuit using sampling.

Tensor network simulators. Qiskit Aer provides a matrix product state (MPS) simulator that can handle certain circuits with up to hundreds of qubits, as long as the entanglement remains limited (low bond dimension). It is not suitable for highly entangled states but works well for circuits with local entanglement structure.

from qiskit_aer import AerSimulator

# Use the matrix_product_state method for larger circuits
sim_mps = AerSimulator(method="matrix_product_state")
# This can handle circuits that would be impossible with statevector simulation,
# provided the entanglement does not grow too large.

Automated Circuit Testing with Assertions

Writing manual debugging checks is fine for one-off investigation, but for circuits you develop iteratively, automated tests catch regressions early. Here is a pattern using Python’s unittest framework.

import unittest
from qiskit import QuantumCircuit
from qiskit.quantum_info import Statevector, DensityMatrix, state_fidelity, entropy, partial_trace
import numpy as np

class TestQuantumCircuits(unittest.TestCase):

    def test_ghz_state_fidelity(self):
        """GHZ state circuit should produce (|000> + |111>)/sqrt(2)."""
        qc = QuantumCircuit(3)
        qc.h(0)
        qc.cx(0, 1)
        qc.cx(1, 2)

        actual = Statevector(qc)
        target = Statevector([1/np.sqrt(2), 0, 0, 0, 0, 0, 0, 1/np.sqrt(2)])
        fid = state_fidelity(actual, target)
        self.assertGreaterEqual(fid, 0.999,
            f"GHZ fidelity too low: {fid:.6f}")

    def test_ghz_entanglement(self):
        """Each qubit in the GHZ state should have entropy log2(2) = 1.0."""
        qc = QuantumCircuit(3)
        qc.h(0)
        qc.cx(0, 1)
        qc.cx(1, 2)

        dm = DensityMatrix(qc)
        for qubit in range(3):
            others = [i for i in range(3) if i != qubit]
            reduced = partial_trace(dm, others)
            ent = entropy(reduced, base=2)
            self.assertGreaterEqual(ent, 0.99,
                f"Qubit {qubit} entropy too low: {ent:.4f}")

    def test_hadamard_maps_zero_to_plus(self):
        """H|0> should equal |+> with fidelity >= 0.999."""
        qc = QuantumCircuit(1)
        qc.h(0)

        actual = Statevector(qc)
        target = Statevector([1/np.sqrt(2), 1/np.sqrt(2)])
        fid = state_fidelity(actual, target)
        self.assertGreaterEqual(fid, 0.999,
            f"H|0> fidelity with |+>: {fid:.6f}")

    def test_hadamard_maps_one_to_minus(self):
        """H|1> should equal |-> with fidelity >= 0.999."""
        qc = QuantumCircuit(1)
        qc.x(0)
        qc.h(0)

        actual = Statevector(qc)
        target = Statevector([1/np.sqrt(2), -1/np.sqrt(2)])
        fid = state_fidelity(actual, target)
        self.assertGreaterEqual(fid, 0.999,
            f"H|1> fidelity with |->: {fid:.6f}")

    def test_parameterized_rotation(self):
        """Ry(theta)|0> should match the analytical result for several angles."""
        for theta in [0, np.pi/6, np.pi/4, np.pi/3, np.pi/2, np.pi]:
            qc = QuantumCircuit(1)
            qc.ry(theta, 0)

            actual = Statevector(qc)
            expected = Statevector([np.cos(theta/2), np.sin(theta/2)])
            fid = state_fidelity(actual, expected)
            self.assertGreaterEqual(fid, 0.999,
                f"Ry({theta:.4f}) fidelity: {fid:.6f}")

# To run these tests, save this code in a file and execute:
# python -m pytest test_circuits.py -v
# or:
# python -m unittest test_circuits.TestQuantumCircuits -v

This testing pattern catches bugs immediately when you modify a circuit. The test_parameterized_rotation test demonstrates how to run the same assertion across multiple parameter values, which helps catch off-by-one errors in rotation angles.

Common Bugs and How to Spot Them

Bug 1: Wrong initial state assumption

Qiskit initializes all qubits to |0>. If you forget to apply an X gate to prepare a qubit in |1>, the circuit runs from the wrong starting point.

# Intended: Bell state starting from |10>
qc_wrong = QuantumCircuit(2)
qc_wrong.h(0)
qc_wrong.cx(0, 1)
# Missing: qc_wrong.x(0) before H

sv_wrong = Statevector(qc_wrong)
print("Wrong initial state amplitudes:", sv_wrong.data)
# Produces |phi+> = (|00> + |11>)/sqrt(2), not (|10> + |01>)/sqrt(2)

Detection: Insert save_statevector(label="initial") as the very first instruction and verify the initial state matches your assumption.

Bug 2: Wrong gate order (non-commuting gates)

# Intended: Ry then Rx
qc_a = QuantumCircuit(1)
qc_a.ry(np.pi / 4, 0)
qc_a.rx(np.pi / 3, 0)

# Accidentally reversed
qc_b = QuantumCircuit(1)
qc_b.rx(np.pi / 3, 0)
qc_b.ry(np.pi / 4, 0)

sv_a = Statevector(qc_a)
sv_b = Statevector(qc_b)
fidelity = state_fidelity(sv_a, sv_b)
print(f"Fidelity between Ry.Rx and Rx.Ry: {fidelity:.6f}")
# Output: Fidelity between Ry.Rx and Rx.Ry: 0.952665
# Will be less than 1.0 because Rx and Ry do not commute

Bug 3: Phase errors invisible in probabilities

Phase errors do not affect individual measurement probabilities but break interference patterns in larger algorithms:

# Correct circuit: H then RZ(pi/4)
qc_correct = QuantumCircuit(1)
qc_correct.h(0)
qc_correct.rz(np.pi / 4, 0)   # T gate phase

# Bugged circuit: accidentally applied RZ(pi/2) instead of RZ(pi/4)
qc_bugged = QuantumCircuit(1)
qc_bugged.h(0)
qc_bugged.rz(np.pi / 2, 0)   # S gate phase (wrong)

sv_c = Statevector(qc_correct)
sv_b = Statevector(qc_bugged)
print(f"Correct state: {sv_c.data}")
print(f"Bugged state:  {sv_b.data}")
print(f"Probabilities match: {np.allclose(sv_c.probabilities(), sv_b.probabilities())}")
print(f"States match:        {np.allclose(sv_c.data, sv_b.data)}")

Probabilities match but states differ. This will cause failures if the circuit is used as a subroutine in a larger algorithm with interference.

Bug 4: Off-by-one in gate parameters

Using pi/4 when you meant pi/2 (or vice versa) is a common typo, especially in rotation gates. The fidelity signature is distinctive: it is close to 1.0 but not exactly 1.0, typically around 0.85 for a pi/4 vs pi/2 error.

# Correct: Ry(pi/2) should rotate |0> to |+>
qc_correct = QuantumCircuit(1)
qc_correct.ry(np.pi / 2, 0)

# Bugged: accidentally wrote pi/4 instead of pi/2
qc_bugged = QuantumCircuit(1)
qc_bugged.ry(np.pi / 4, 0)

sv_correct = Statevector(qc_correct)
sv_bugged = Statevector(qc_bugged)

fid = state_fidelity(sv_correct, sv_bugged)
print(f"Fidelity: {fid:.6f}")
# Output: approximately 0.853553
# This "close but not perfect" fidelity is the signature of a parameter error.

# Compare amplitudes to see the exact difference
print(f"Correct: {sv_correct.data}")
# Output: [0.70710678, 0.70710678]
print(f"Bugged:  {sv_bugged.data}")
# Output: [0.92387953, 0.38268343]

Detection signature: Fidelity between 0.7 and 0.99 (not near 0 and not near 1) usually indicates a parameter error. The state is “in the right neighborhood” but not exactly correct.

Bug 5: Wrong qubit index in two-qubit gate

Swapping the control and target of a CNOT gate, or applying it to the wrong qubit pair, produces an entanglement pattern in the wrong basis.

# Correct: CNOT from qubit 0 (control) to qubit 1 (target)
qc_correct = QuantumCircuit(3)
qc_correct.h(0)
qc_correct.cx(0, 1)  # Entangle qubits 0 and 1

# Bugged: CNOT applied to wrong qubit pair
qc_bugged = QuantumCircuit(3)
qc_bugged.h(0)
qc_bugged.cx(0, 2)  # Accidentally entangled qubits 0 and 2

sv_correct = Statevector(qc_correct)
sv_bugged = Statevector(qc_bugged)

print("Correct circuit amplitudes:")
for i, amp in enumerate(sv_correct):
    if abs(amp) > 1e-10:
        print(f"  |{format(i, '03b')}> : {amp:.4f}")
# Output: |000> and |110> have amplitude 0.7071 (qubits 0 and 1 entangled)

print("\nBugged circuit amplitudes:")
for i, amp in enumerate(sv_bugged):
    if abs(amp) > 1e-10:
        print(f"  |{format(i, '03b')}> : {amp:.4f}")
# Output: |000> and |101> have amplitude 0.7071 (qubits 0 and 2 entangled)

Detection signature: The nonzero amplitudes appear in different basis states than expected. Looking at which qubits are correlated in the output reveals which qubit index was wrong.

Bug 6: Missing gate

Forgetting to apply a gate is easy to do, especially in longer circuits. The save_statevector technique at multiple points isolates exactly where the state diverges from the expected evolution.

from qiskit_aer import AerSimulator

# Correct: H, CNOT, then Z on qubit 0
qc_correct = QuantumCircuit(2)
qc_correct.h(0)
qc_correct.cx(0, 1)
qc_correct.z(0)

# Bugged: forgot the Z gate
qc_bugged = QuantumCircuit(2)
qc_bugged.h(0)
qc_bugged.cx(0, 1)
# Missing: qc_bugged.z(0)

# Add probes to the bugged circuit to find where state diverges
qc_probed = QuantumCircuit(2)
qc_probed.h(0)
qc_probed.save_statevector(label="step1")
qc_probed.cx(0, 1)
qc_probed.save_statevector(label="step2")
# The missing Z gate means there is no step3

sim = AerSimulator()
result = sim.run(qc_probed).result()

# Compare each probe against expected state at that point
expected_step1 = Statevector(QuantumCircuit(2).compose(
    QuantumCircuit(2, name="h").compose(
        QuantumCircuit(1).compose(QuantumCircuit(1)), qubits=[0]
    )
))

# Simpler approach: build expected states directly
qc_exp1 = QuantumCircuit(2)
qc_exp1.h(0)
expected_step1 = Statevector(qc_exp1)

qc_exp2 = QuantumCircuit(2)
qc_exp2.h(0)
qc_exp2.cx(0, 1)
expected_step2 = Statevector(qc_exp2)

qc_exp3 = QuantumCircuit(2)
qc_exp3.h(0)
qc_exp3.cx(0, 1)
qc_exp3.z(0)
expected_step3 = Statevector(qc_exp3)

sv1 = Statevector(result.data(0)["step1"])
sv2 = Statevector(result.data(0)["step2"])

print(f"Step 1 fidelity: {state_fidelity(sv1, expected_step1):.6f}")
# Output: 1.000000 (H gate correct)
print(f"Step 2 fidelity: {state_fidelity(sv2, expected_step2):.6f}")
# Output: 1.000000 (CNOT correct)
print(f"Final fidelity:  {state_fidelity(sv2, expected_step3):.6f}")
# Output: less than 1.0 (state matches step 2 but not the expected final state)
# This tells you the bug is AFTER step 2: the Z gate is missing.

Detection signature: The state matches all intermediate expected states up to a certain point, then diverges. The divergence point is where the missing gate should have been applied.

Bug 7: Applying a gate twice (double application)

Some gates are their own inverse. The Hadamard gate satisfies HH = I, so applying it twice cancels the effect entirely. The X gate similarly satisfies XX = I. If you accidentally apply one of these gates twice, the circuit behaves as if the gate was never applied.

# Correct: apply H once to get |+>
qc_correct = QuantumCircuit(1)
qc_correct.h(0)

# Bugged: accidentally applied H twice (copy-paste error)
qc_bugged = QuantumCircuit(1)
qc_bugged.h(0)
qc_bugged.h(0)  # This cancels the first H

sv_correct = Statevector(qc_correct)
sv_bugged = Statevector(qc_bugged)

print(f"Correct (single H): {sv_correct.data}")
# Output: [0.70710678, 0.70710678]  (|+> state)
print(f"Bugged (double H):  {sv_bugged.data}")
# Output: [1.+0.j, 0.+0.j]  (back to |0>, as if H was never applied)

fid = state_fidelity(sv_correct, sv_bugged)
print(f"Fidelity: {fid:.6f}")
# Output: 0.500000

Detection signature: A fidelity of exactly 0.5 between the actual and expected states often indicates a self-inverse gate applied twice (or not at all). The state looks like the gate was never applied, because the gate canceled itself out.

This bug is especially common with copy-paste errors in circuit construction, or when refactoring a circuit and accidentally leaving a gate that should have been removed.

Key Takeaways

  • Statevector(qc) gives you immediate access to all amplitudes without running a simulator job.
  • save_statevector(label=...) is the quantum equivalent of a print-statement debugger.
  • state_fidelity() provides a single number to assert circuit correctness and correctly ignores global phase.
  • DensityMatrix handles noisy and mixed states where Statevector cannot.
  • partial_trace and entropy let you inspect individual qubits and detect entanglement in multi-qubit systems.
  • Qiskit uses little-endian qubit ordering: qubit 0 is the rightmost bit. Write a helper function to avoid confusion.
  • Phase errors are invisible in measurement probabilities but detectable by comparing statevectors directly.
  • Statevector debugging scales to approximately 25 qubits on consumer hardware. Beyond that, use sampling-based tests or tensor network simulators.
  • Automated tests using state_fidelity assertions catch regressions early and document expected circuit behavior.
  • Common bug signatures: fidelity near 0.85 suggests a parameter error, fidelity of 0.5 suggests a doubled self-inverse gate, and fidelity near 0 suggests an orthogonal state (wrong phase or wrong gate entirely).

Was this tutorial helpful?