How Do Tensor Networks Simulate Quantum Circuits?

Tensor network contraction algorithms enable classical supercomputers to simulate shallow or noisy quantum circuits by representing quantum states and gates as interconnected multidimensional arrays rather than a single exponentially large state vector. Instead of tracking every possible quantum state simultaneously, these algorithms exploit the limited entanglement present in shallow or noisy circuits to compress the computation and contract the network in an optimized order.

The Exponential Bottleneck of Classical Simulation

A quantum system of \(n\) qubits is described by a state vector in a \(2^n\)-dimensional Hilbert space. Storing this vector directly on a classical computer requires memory that scales as \(O(2^n)\). For \(n = 50\) qubits, even storing a single state vector with double-precision complex numbers demands over \(16\) petabytes of memory, quickly exhausting the capacity of modern supercomputers.

However, most quantum algorithms—especially those run on near-term noisy intermediate-scale quantum (NISQ) hardware—do not explore the entire Hilbert space uniformly. When a circuit is shallow (having few layers of gates) or noisy (suffering from decoherence and gate errors), the amount of quantum entanglement between distant qubits remains bounded. Tensor network methods take advantage of this structure.

Representing Quantum Circuits as Tensor Networks

In a tensor network, every component of a quantum circuit is mapped to a tensor—a generalization of scalars, vectors, and matrices to higher dimensions:

  • Qubit initializations (such as \(|0\rangle\)) are rank-1 tensors (vectors).
  • Single-qubit gates are rank-2 tensors (matrices) with one input index and one output index.
  • Two-qubit gates (such as CNOT or CZ) are rank-4 tensors with two input indices and two output indices.
  • Measurements are rank-1 tensors that project onto a basis state.

When gates are applied sequentially to qubits, the corresponding tensors are connected along shared indices, forming a graph. Computing a specific measurement probability or amplitude \(\langle x | U | 0 \rangle^{\otimes n}\) corresponds to contracting all internal indices of this tensor network—summing over all shared indices until only a single scalar remains.

Why Contraction Order Matters

The computational cost of contracting a tensor network depends heavily on the order in which tensors are multiplied together. If two tensors share an index of dimension \(d\) (for qubits, \(d = 2\)), contracting them eliminates that index and produces a new tensor whose rank is the sum of the remaining uncontracted indices.

Finding the optimal contraction order is equivalent to finding a contraction tree with the minimum treewidth (or minimum maximum intermediate tensor size). A poor contraction order can cause intermediate tensors to grow to size \(2^n\), recovering the exponential bottleneck. An optimized contraction order, by contrast, keeps intermediate tensors small by contracting along directions with low entanglement—such as contracting across time (circuit depth) rather than space (number of qubits) when the circuit is shallow.

Supercomputers use heuristic graph-partitioning and dynamic programming algorithms to discover near-optimal contraction paths before executing the actual floating-point operations.

Exploiting Shallow Depth and Noise

Tensor network contraction excels in two regimes common to current quantum hardware:

  1. Shallow circuits: If a circuit has depth \(D\) on a 2D grid of qubits, the entanglement across any cut grows at most proportionally to \(D\). By contracting the network along the spatial dimensions first or slicing over a small number of indices, the memory footprint scales roughly as \(O(2^{c \cdot D})\) for a small constant \(c\), rather than \(O(2^n)\).
  2. Noisy circuits: Quantum noise (such as depolarizing or amplitude-damping channels) degrades long-range entanglement. In tensor network representations like Matrix Product States (MPS) or Projected Entangled Pair States (PEPS), noise allows algorithms to truncate small singular values during Singular Value Decomposition (SVD). Dropping singular values below a threshold \(\epsilon\) introduces a controlled approximation error while drastically reducing the bond dimension \(\chi\) needed to represent the state.

Scaling Across Supercomputers

Once an efficient contraction path is found, the workload maps naturally onto distributed supercomputing architectures:

  • Index slicing: To fit intermediate tensors into the memory of individual GPUs or compute nodes, large indices are "sliced" (split into independent sums). Each slice can be contracted independently in parallel without inter-node communication, and the final results are summed at the end.
  • Tensor-core acceleration: Modern supercomputers equipped with GPUs perform tensor contractions as batched matrix multiplications (GEMM), achieving high utilization of specialized hardware like Tensor Cores.

By combining optimal contraction ordering, index slicing, and entanglement truncation, tensor network algorithms allow classical supercomputers to simulate quantum circuits with dozens or even hundreds of qubits—provided the circuits remain sufficiently shallow or noisy.