7. Numerical Computation
7.1Floating-point arithmetic
Floating-point numbers are scientific notation in binary: a sign, a mantissa, and an exponent, giving about 15–16 significant decimal digits for float64. They are not real numbers: most values (0.1 among them) have no exact representation, arithmetic rounds at every step, and familiar algebra partially fails — addition is not associative, (a + b) + c ≠ a + (b + c) in general. The constants of the craft: machine epsilon (≈ 2.2×10⁻¹⁶ for float64, np.finfo(float).eps) is the relative granularity; comparisons use tolerances (np.allclose, np.isclose — never ==); catastrophic cancellation amplifies relative error when subtracting nearly equal numbers (compute |⟨φ|ψ⟩|² as the squared modulus, not as ⟨φ|ψ⟩·⟨φ|ψ⟩ with a separately-computed phase). Quantum states compound this by the million: 2ⁿ amplitudes each rounding per gate. Nothing here is exotic — it is the same arithmetic every spreadsheet runs — but quantum simulation multiplies its failures by 2ⁿ.
7.2Numerical error
Every numerical result carries two error layers: roundoff (representation and arithmetic granularity) and truncation (approximation built into an algorithm — a truncated series, an iterative solve stopped early). Engineering means budgeting both. Roundoff accumulates gradually: a circuit of 10⁴ gates on float64 amplitudes accumulates relative error around 10⁴·10⁻¹⁶ ≈ 10⁻¹² — negligible per gate, visible in aggregate, and the reason simulators re-normalize states periodically rather than trusting the algebra to keep ‖ψ‖ = 1. Truncation error is algorithmic and usually dominant: a Krylov exponential stopped at m steps, a Trotter split of order Δt², a finite-difference derivative with step h. The discipline: for every approximation, know its order — how error scales with the knob (steps, threshold, h) — and verify empirically by halving the knob and watching the error fall by the predicted factor. An error estimate you have not validated by convergence testing is decoration.
7.3Conditioning
Conditioning measures how violently a problem's answer responds to small input perturbations — a property of the problem, not the algorithm. The condition number κ(A) = σ_max/σ_min (ratio of largest to smallest singular value) prices linear solves: relative error in x of Ax = b can reach κ times the relative error in the inputs, so κ = 10¹⁰ on float64 data means you may retain only ~6 correct digits no matter how good the solver. Check it: np.linalg.cond(A); act when it warns — rescale variables, regularize, or reformulate. Quantum instances: density matrices near pure states are rank-deficient and ill-conditioned (regularize before inverting covariance matrices in tomography); eigenvectors of nearly degenerate eigenvalues are individually unstable even though the invariant subspace is fine (compare subspaces, not vectors, 4.16); and least-squares fits of noisy calibration data need ridge regularization or they chase noise. Rule: ill-conditioning is a property to be designed around, never merely endured.
7.4Numerical linear algebra
Numerical linear algebra is the discipline of computing matrix answers stably and fast, and it is the load-bearing wall of every quantum simulator. The core decompositions, each with a specialty: LU (with pivoting) for general solves — np.linalg.solve; Cholesky for symmetric/Hermitian positive-definite systems, cheaper and stable; QR for orthonormalization and least squares — also how simulators re-orthonormalize state batches (3.19); SVD, the most informative, giving singular values, ranks, conditioning, and best low-rank approximation in one call — np.linalg.svd; and eigendecompositions eigh/eig for spectra (4.17). Cost is O(n³) for dense n×n across the board, which sets the practical ceiling: dense state-vector simulation tops out near n ≈ 25–28 qubits on a laptop, n ≈ 40 on a large machine. The engineering habit: pick the decomposition matching your matrix's structure (Hermitian? positive-definite? sparse?) — the wrong generic routine costs accuracy and an order of magnitude in time.
7.5Sparse matrices
A matrix is sparse when most entries are zero, and quantum operators are sparsity royalty: an n-qubit Hamiltonian with local interactions has O(n·polylog) nonzero entries in a 2ⁿ×2ⁿ matrix — storing it densely is malpractice past ~15 qubits. scipy.sparse (csr_matrix, csc_matrix) stores (row, column, value) triples and computes with them: sparse matvec costs O(nnz), not O(n²), and applying I⊗…⊗G⊗…⊗I to a state vector touches only the entries the gate couples — the trick that pushes laptop simulation to 30+ qubits. Sparse-aware eigensolvers (scipy.sparse.linalg.eigsh, Lanczos-based) extract a few extremal eigenvalues of huge matrices without forming them. Two caveats: sparsity is per-operator and degrades under multiplication — products of sparse matrices fill in, so Trotter steps apply factors sequentially instead of composing them; and indices remain 64-bit, so 2ⁿ exceeding 2⁶³ is a wall even sparse code cannot pass. Structure beats size; exploit both.
7.6Eigenvalue algorithms
Production eigenvalue computation is iterative: dense solvers (LAPACK behind numpy's eigh) reduce to tridiagonal form then iterate (QR algorithm) — robust, O(n³), all eigenvalues at once. But quantum problems usually want a few eigenvalues of a huge matrix, and for that the iterative family rules. Lanczos (Hermitian case) and Arnoldi (general) build a small Krylov subspace from repeated matrix-vector products — needing only the ability to apply A, never to store it — and extract Ritz-value approximations of extremal eigenpairs; scipy.sparse.linalg.eigsh(A, k=6, which='SA') is the workhorse. Convergence is fastest for well-separated extremal eigenvalues and slow inside spectra clusters; preconditioning and shift-invert strategies (which do need solves, hence good conditioning) fix the hard cases. The same Krylov machinery computes e^{A}v for large sparse A — matrix exponentials of Hamiltonians without diagonalization, the standard method for quantum dynamics at sizes where eigh is impossible.
7.7Optimization
Optimization — minimizing a scalar function f(θ) over parameters — powers the classical half of variational quantum algorithms (chapter 46): the parameters θ of a parameterized circuit are tuned by a classical optimizer that only sees sampled costs. The landscape of methods: gradient descent and its momentum/adaptive variants (Adam) when gradients exist — on quantum hardware, gradients come from parameter-shift rules, exact but costing two circuit evaluations per parameter; derivative-free methods (Nelder–Mead, COBYLA) when noise makes gradients meaningless; and global/heuristic methods (SPSA — perturb-and-difference, famously robust to noisy cost evaluations) for barren, flat, or deceptive landscapes. Every method carries tuning folklore: learning rates, initial points, stopping criteria — and every noisy quantum cost function adds variance on top (6.5), so optimizer convergence must itself be assessed statistically. The unifying discipline from calculus: check gradients numerically (scipy.optimize.check_grad-style finite differences) before trusting any optimizer's complaints about them.
7.8Automatic differentiation
Automatic differentiation (AD) computes exact derivatives of code: not symbolic expansion, not finite-difference approximation, but the chain rule applied mechanically to every elementary operation, giving gradients to machine precision at a small constant factor of runtime. Two modes: reverse mode (backpropagation) computes all partial derivatives of one output in one sweep — cost roughly independent of parameter count, hence deep learning's engine; forward mode is cheaper for few inputs, many outputs. The quantum connection is direct: autodiff simulators differentiate entire circuits, making hybrid quantum-classical training gradient-based end to end; and the parameter-shift rule (7.7) is AD's hardware analogue — for gates of the form e^{iθP}, ∂f/∂θ = [f(θ + π/2) − f(θ − π/2)]/2 exactly, sampled on real devices. Tooling: JAX and PyTorch for classical graphs, qiskit's and PennyLane's built-in AD for hybrid workflows. The engineer's rule: if your loss is differentiable, hand-rolled finite differences are waste — and usually wrong.
7.9Monte Carlo simulation
Monte Carlo reappears here as a numerical method with engineering teeth: when a quantity has no closed form — an integral over configurations, a noise-averaged observable, a tail probability — simulate the underlying randomness N times and average. The machinery is chapter 6's (sampling, √N convergence, error bars), applied to physics: simulate noisy circuits shot-by-shot with stochastic gate errors and decoherence events; estimate detection probabilities of rare error syndromes by importance sampling that forces syndrome events and re-weights; average observable estimates over draws of unknown noise parameters. Variance reduction is the craft: importance sampling (sample the rare, reweight), control variates (subtract a correlated quantity with known mean), and antithetic draws (pair each sample with its mirror) routinely buy one to three orders of magnitude in shots. In error-correction research (chapter 37), Monte Carlo over syndrome histories is the standard evaluation of decoders — and it is embarrassingly parallel, so your laptop's cores and a concurrent.futures map go a long way.
7.10Performance considerations
Simulation performance is mostly memory bandwidth and algorithm choice, rarely raw FLOPS. The hierarchy of wins, in order of leverage: algorithmic (avoid forming 2ⁿ×2ⁿ operators — apply gates factor-wise, 7.11; exploit sparsity, 7.5; use stabilizer or tensor-network representations when the circuit permits, 5.14); precision (float32 halves memory and doubles throughput versus float64, tolerable for sampling experiments though not for phase-sensitive accumulation — know which your experiment is); vectorization (numpy applies gates to whole state batches per call — a Python loop over amplitudes is a 100× self-inflicted tax); parallelism (multicore via vectorized BLAS threads or multiprocessing, since state updates parallelize trivially across amplitudes); and memory layout (contiguous complex128 arrays, minimal copies — np.kron chains allocate massively; build operators in-place where you can). Profile before optimizing — cProfile plus a stopwatch on the hot gate-application loop — and expect the answer to be "you materialized something you should have applied on the fly".
7.11Classical simulation of quantum systems
This is the chapter's summit: the techniques above assemble into a working n-qubit simulator, the tool you will use more than any other in this book. The dense state-vector engine is 200 lines: states as complex128 arrays of length 2ⁿ; gates as small matrices applied via reshape-and-broadcast (apply single-qubit G on qubit k by reshaping the state to shape (2, 2, …, 2) and tensordoting along axis k — no np.kron, no 2ⁿ×2ⁿ operator ever formed); measurements by sampling from |amplitudes|²; noise by density matrices (memory ×2², halving your qubit budget) or by stochastic trajectories (sample which error occurred, evolve the state vector — Monte Carlo over noise realizations, 7.9). Beyond dense: stabilizer simulators run certain circuits (Clifford gates, Pauli measurements) in polynomial time and simulate error-correction experiments at hundreds of qubits; tensor networks buy entanglement-limited sizes. Chapter 14 builds the engine properly; this section is its numerical foundation.
7.12Why simulation becomes exponentially difficult
The closing honesty: no cleverness removes the exponential, because the exponential is the information. A generic n-qubit state has 2ⁿ complex amplitudes of independent content; any representation that stores fewer must be discarding structure, and works exactly where that structure is absent or shallow. The escape hatches are all structure-conditional: product and low-entanglement states compress (5.7, 5.14); Clifford circuits simulate in polynomial time but are efficiently classically simulable precisely because they cannot do universal computation; sparse operators tame gates but not generic states. Random deep circuits — the very ones used for supremacy claims — are designed to evade all of these, which is why simulating 50+ qubits of them strains the world's best machines and why each quantum announcement is followed by a classical counterattack with better networks. For your practice: know which regime a problem sits in before choosing a tool, and treat "simulate it" as a research decision, not a default.