17. Building a Quantum Simulator

17.1Why Build One Yourself

You will use Aer for everything real, so why write 300 lines of simulator? Because a simulator is the one quantum system you can see all the way into — every amplitude, every phase — and building one converts Part V's mathematics from notation into muscle memory. The exercise is the field's equivalent of writing an interpreter: after it, "the simulator applies the gate" is no longer a black box, and Aer's method flags (16.6) become engineering choices you understand. Budget: one focused weekend. Requirements: numpy, nothing else. The project spec at this chapter's end defines acceptance tests — write them first.

17.2State-Vector Simulation

The core loop: maintain the state as a complex numpy array of shape (2ⁿ,), apply each gate by updating amplitudes, and at the end sample from |amplitudes|². For n qubits, index i = Σ bₖ·2ᵏ encodes the computational basis state |bₙ₋₁…b₀⟩. A gate on qubit k mixes only pairs of amplitudes differing in bit k — the observation that makes simulation O(2ⁿ) per gate instead of O(4ⁿ). Memory is the binding constraint: 128 bytes per amplitude (complex128) means 30 qubits = 137 GB. Everything else in this chapter elaborates this loop.

17.3Single-Qubit Simulation

import numpy as np

I2 = np.eye(2, dtype=complex)
H = np.array([[1, 1], [1, -1]], dtype=complex) / np.sqrt(2)
X = np.array([[0, 1], [1, 0]], dtype=complex)

def apply_1q(state, gate, k, n):
    """Apply `gate` to qubit k of an n-qubit state vector."""
    gk = np.kron(np.kron(np.eye(2**k), gate), np.eye(2**(n-k-1)))
    return gk @ state

This is correct and pedagogically honest — and hopelessly slow (O(4ⁿ) per gate). Run it for n=2..20, time it, and keep the numbers; 17.7 will beat them by orders of magnitude, and you should understand exactly why the first version dies: it multiplies by a mostly-identity 2ⁿ×2ⁿ matrix.

17.4Matrix Multiplication

The naive approach treats the whole circuit as one matrix: U_total = U_m···U_2·U_1, then U_total @ state. Full 2ⁿ×2ⁿ matrices are infeasible past n≈13 (a 27 GB dense matrix), so serious simulators never build them. But the composition view is worth keeping for n≤10 as a ground truth: build U_total with np.kron and single products, apply once, and use the result to unit-test the fast implementation (17.7). Differential testing against the naive version is the single most effective bug-catcher in this project — keep both implementations in the test suite forever.

17.5Multi-Qubit Simulation

Controlled gates are the first real challenge. CX on (control c, target t): for each pair of amplitudes (i, i XOR 2ᵗ) where bit c of i is 1, swap-mix them: (a_i, a_j) → (a_j, a_i) for X; more generally apply the 2×2 gate to pairs conditioned on the control bit. The mask idiom: mask_c = 1 << c; for i in range(2**n): if i & mask_c: pair with j = i ^ (1 << t). Watch for double-counting — iterate over pairs, not all indices. This loop, vectorized correctly, is 90% of the project.

17.6Tensor Products

np.kron is how operators on subsystems become operators on the whole: H on qubit 0 of a 3-qubit system is kron(H, I, I) — with the convention question: does qubit 0 get the leftmost factor? Decide once (numpy row-major favors leftmost = most significant bit; Qiskit's little-endian favors the opposite) and write a test that documents your choice, because a transposed endianness bug produces plausible-looking wrong answers. Tensor products also give you product-state initialization and the marginalization needed for measurement simulation. Part V's math becomes concrete here.

17.7Applying Gates Efficiently

The production-grade approach: reshape + einsum. View the state as an n-dimensional tensor of shape (2,)*n with axis k = qubit k (choose your convention), then apply a 2×2 gate to axis k via np.tensordot(gate, state, axes=([1],[k])) and move the axis back with np.moveaxis. Cost: O(2ⁿ) per gate — optimal — with numpy doing the inner loops in C. Two-qubit gates: tensordot with a 4×4 matrix over two axes, or decompose into O(1) single-axis ops. Benchmark: your n=26 simulation should apply gates in well under a second; if not, you are still materializing 2ⁿ×2ⁿ matrices somewhere.

17.8Measurement Simulation

Given final amplitudes, probabilities are np.abs(state)2 (your Born-rule implementation, already built in Chapter 1 — reuse it). Sampling one shot: np.random.choice(2n, p=probs); N shots: vectorize or loop. Post-measurement state (needed for mid-circuit measurement): zero out amplitudes inconsistent with the outcome and renormalize. Decide your convention for mapping the sampled integer to a bit-string (and test it!). Measurement is also where memory pressure spikes: precompute probabilities once, sample many times, and never materialize the full 2ⁿ probability array more than once per measurement point.

17.9Random Sampling

Simulators owe you two kinds of randomness, done right. Outcome sampling: use a seeded np.random.Generator (PCG64) for reproducibility — your property tests (16.15) depend on it. And honest noise: a plain state-vector simulator is ideal; adding realistic noise means Kraus-channel sampling (Part IX will build this on top of your engine — design the API for it now: apply_channel(state, kraus_ops) -> state). Also implement exact expectation values ⟨ψ|O|ψ⟩ via np.vdot(state, O @ state) for the cases where sampling is wasteful. A simulator that only samples is half a simulator.

17.10Circuit Parsing

Give your engine a front end so tests read like circuits. Define a tiny text format:

H 0
CX 0 1
T 1
M 0 1

Parse lines → a list of (gate, operands) instructions; validate (qubit indices in range, no measurement mid-circuit unless you support it); then interpret. Round-trip test: convert a Qiskit circuit to this format (qc.draw('text') parsing is overkill — use qiskit.qasm2 or export your own), run on both engines, assert equal distributions with the same seed. Now your simulator is a dual of Qiskit's, and every future chapter can check one against the other.

17.11Performance Optimization

Ordered by payoff. One: tensordot over full-matrix apply (100–1000×). Two: dtype=complex128 consistently — silent upcasting to complex256 or object arrays is a classic slowdown. Three: fuse single-qubit gate sequences per qubit before applying (matrix product of 2×2s, one tensordot). Four: batch shots by sampling once per distribution. Five: np.einsum with optimize=True for multi-axis contractions. Profile before optimizing (cProfile, line_profiler) — in simulation code the profile is almost never where intuition points. Realistic target: ~10⁶ gates/second at n=25 on a laptop M-series chip; Aer will still beat you 10×, and that's fine.

17.12Memory Complexity

The state vector is 16·2ⁿ bytes; the wall arrives at n=30 (172 GB) on commodity hardware, n≈34 on a 512 GB machine. Be precise about what "simulating n qubits" costs: time O(gates·2ⁿ), memory O(2ⁿ) — and notice the asymmetry with hardware, which costs O(n) qubits to hold the same state. Extend Chapter 1's memory experiment: your own simulator, n = 24…30, plot time-per-gate (flat) and max memory (exponential). The exponential you measure on your laptop is the same constant that protects a 50-qubit state from simulation by the biggest supercomputer — the field's honest boundary, experienced firsthand.

17.13Sparse Simulation

Many important states are sparse: computational basis states, and states reachable by circuits with few non-Clifford gates. Sparse simulators store (index, amplitude) pairs in a dict or hash map, and a gate on qubit k touches at most 2ⁿ⁻¹ occupied pairs. When does this win? Controlled additions (Shor-style arithmetic) create states with O(poly) occupied basis states per intermediate step — sparse simulation of Shor N=15 state-preparation runs in kilobytes. Implement SparseSimulator sharing your parser (17.10), and test it against the dense engine on small n. Qiskit's Aer doesn't ship one; research codes (and Google's 2019 Supremacy crossover claims) live in this regime.

17.14Stabilizer Simulation

The Gottesman–Knill theorem: circuits made only of {H, S, CX, Paulis, measurement of Z} are classically simulable in O(n²) time and memory — no state vector at all. The state is represented by n stabilizer generators (Pauli strings with signs); each Clifford gate updates the tableau in polynomial time. This is why thousands of qubits of error-correcting circuitry (Part X: surface codes are Clifford circuits!) simulate fine while 50 algorithmic qubits don't. Implement a mini-tableau (Chase Gordon's stim is the reference — install it, read its paper) and verify: your stabilizer engine agrees with the dense engine on Clifford circuits of 10–30 qubits, in microseconds vs. seconds.

17.15Tensor-Network Simulation

The frontier of classical simulation: represent the state as a network of small tensors (matrix product states, PEPS) whose size tracks entanglement, not qubit count. Low-entanglement circuits (1D dynamics, shallow depths) simulate to 100+ qubits; highly entangled circuits hit the bond-dimension wall — which is exactly how Google argued its 53-qubit 2019 result was hard to simulate classically (and how IBM's counter-analysis argued the crossover was reachable). You won't implement MPS here (Part XVIII offers it as a project), but know quimb and ITensor, and know the rule: tensor networks trade memory for entanglement assumption — the perfect counterpart to 17.12's hard exponential.