Case Study: Building a State-Vector Simulator from Tensor Products

Executive Summary

The fastest way to stop being confused by tensor products is to write a simulator. In roughly forty lines of NumPy you can build something that runs real quantum circuits, and in doing so you will confront every piece of Chapter 3's machinery as an engineering constraint rather than a definition: why the state vector has $2^n$ entries, why a single-qubit gate on qubit $k$ of an $n$-qubit register is not a $2\times2$ matrix, why the naive implementation dies at 20 qubits, and how to fix it.

Skills applied

• Constructing $n$-qubit states and operators as tensor products (§3.11).
• Embedding a local operator into a global Hilbert space via $I^{\otimes a}\otimes U \otimes I^{\otimes b}$ (§3.12).
• Applying operators by reshaping rather than by matrix multiplication.
• Reasoning about the exponential memory wall quantitatively.

Phase 1: The state

An $n$-qubit pure state is a complex vector of length $2^n$, indexed by the integer whose binary expansion is the basis-state label:

import numpy as np

def zero_state(n):
    psi = np.zeros(2**n, dtype=complex)
    psi[0] = 1.0                     # |00...0>
    return psi

Index 6 in a 3-qubit register is $|110\rangle$. That is the entire encoding — the tensor product is implicit in the indexing, which is the first genuinely useful realization.

Phase 2: The naive gate application

To apply a single-qubit gate $U$ to qubit $k$ of $n$, the textbook construction is

$$U_k = I^{\otimes k} \otimes U \otimes I^{\otimes (n-k-1)}$$

def embed(U, k, n):
    op = np.array([[1]], dtype=complex)
    for q in range(n):
        op = np.kron(op, U if q == k else np.eye(2))
    return op

def apply_naive(psi, U, k, n):
    return embed(U, k, n) @ psi

This is correct and it is how the mathematics is written. It is also unusable.

The cost

For $n$ qubits, embed builds a $2^n \times 2^n$ matrix — $4^n$ complex numbers at 16 bytes each:

The state is affordable well past 20 qubits. The operator is not affordable past about 14. We are storing an enormous matrix that is almost entirely identity structure.

Finding. The tensor product is the right mathematics and the wrong data structure. This is the single most important practical lesson of Chapter 3.

Phase 3: The reshape trick

A single-qubit gate touches one index. Reshape the state vector so that index is exposed, contract over it, and reshape back — never materializing the big matrix:

def apply(psi, U, k, n):
    psi = psi.reshape([2] * n)                  # tensor of rank n
    psi = np.tensordot(U, psi, axes=([1], [k])) # contract U's input with qubit k
    psi = np.moveaxis(psi, 0, k)                # restore axis order
    return psi.reshape(2**n)

Cost: $O(2^n)$ time, $O(2^n)$ memory. No $4^n$ anywhere. The same trick generalizes to two-qubit gates by contracting over two axes.

def apply2(psi, U, k, l, n):
    psi = psi.reshape([2] * n)
    U = U.reshape(2, 2, 2, 2)                   # (out_k, out_l, in_k, in_l)
    psi = np.tensordot(U, psi, axes=([2, 3], [k, l]))
    psi = np.moveaxis(psi, [0, 1], [k, l])
    return psi.reshape(2**n)

This is, in essence, how every serious state-vector simulator works.

Phase 4: Measurement

Born rule, directly:

def probabilities(psi):
    return np.abs(psi)**2

def sample(psi, shots=1024, rng=None):
    rng = rng or np.random.default_rng()
    n = int(np.log2(len(psi)))
    outcomes = rng.choice(len(psi), size=shots, p=probabilities(psi))
    counts = {}
    for o in outcomes:
        key = format(o, f'0{n}b')
        counts[key] = counts.get(key, 0) + 1
    return counts

Note what this makes obvious: the simulator has the whole amplitude vector in hand and still must sample to produce output, because that is what measurement does. Having the state is not having the answer.

Phase 5: Running a real circuit

Bell state, then a three-qubit GHZ:

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

psi = zero_state(2)
psi = apply(psi, H, 0, 2)
psi = apply2(psi, CNOT, 0, 1, 2)
print(np.round(psi, 3))     # [0.707+0j, 0, 0, 0.707+0j]  ->  (|00> + |11>)/sqrt(2)

psi = zero_state(3)
psi = apply(psi, H, 0, 3)
psi = apply2(psi, CNOT, 0, 1, 3)
psi = apply2(psi, CNOT, 1, 2, 3)
print(sample(psi, 1000))    # ~{'000': 500, '111': 500}

The GHZ output is worth staring at. Only $|000\rangle$ and $|111\rangle$ appear, each about half the time. No amount of examining qubit 0 alone predicts qubit 2 — yet they always agree. That correlation is the thing the tensor-product formalism exists to describe, and it is why the state does not factor.

Phase 6: Verifying non-factorizability

The formal test for whether a two-qubit state is a product state is the rank of its coefficient matrix. Write $|\psi\rangle = \sum_{ij} c_{ij}|ij\rangle$ and arrange $c$ as a $2\times2$ matrix:

def is_product_state(psi2, tol=1e-10):
    C = psi2.reshape(2, 2)
    return np.linalg.matrix_rank(C, tol=tol) == 1

Rank 1 means $C = ab^T$ for vectors $a, b$, i.e. $|\psi\rangle = |a\rangle\otimes|b\rangle$. For the Bell state, $C = \frac{1}{\sqrt2}\begin{pmatrix}1&0\\0&1\end{pmatrix}$ — rank 2, not a product. The singular values of $C$ are the Schmidt coefficients, and the number of nonzero ones is the Schmidt rank, the cleanest single measure of how entangled a bipartite pure state is.

Discussion Questions

1. The naive implementation is a faithful transcription of the mathematics and is unusable by 15 qubits. What does that suggest about reading formalism as implementation guidance generally?
2. The reshape trick costs $O(2^n)$ per gate. For a circuit of $g$ gates the total is $O(g2^n)$. At what qubit count does a 1,000-gate circuit stop fitting in an hour on a laptop?
3. Why can a simulator holding the exact state vector still not tell you the measurement outcome in advance?
4. Tensor-network simulators beat state-vector simulators on some circuits by orders of magnitude. What property of a circuit would you guess makes it tensor-network-friendly?

Your Turn: Extensions

• Add $R_x, R_y, R_z$ and a controlled-$U$ constructor; run the Chapter 1 Grover example on 4 qubits.
• Implement mid-circuit measurement: project onto the outcome, renormalize, and continue.
• Benchmark apply_naive against apply for $n = 8, 10, 12, 14$ and plot the crossover.
• Compute the Schmidt coefficients of a randomly generated two-qubit state and confirm they are the singular values of the reshaped coefficient matrix.

Key Takeaways

• An $n$-qubit state is a length-$2^n$ vector whose index is the basis-state bit string; the tensor product lives in the indexing.
• Never materialize $I^{\otimes a}\otimes U\otimes I^{\otimes b}$: reshape and contract instead. $O(2^n)$ beats $O(4^n)$ decisively.
• The memory wall for state-vector simulation is the state itself: roughly 30–40 qubits on the largest machines, permanently.
• Entanglement is non-factorizability, and it is detectable as the rank of the reshaped coefficient matrix.
• Holding the full quantum state still leaves you sampling for answers — simulation reproduces measurement, it does not evade it.