4. Matrices
4.1Matrix representation
A matrix is a rectangular array of numbers with defined operations; a size-m×n matrix has m rows and n columns and maps n-dimensional vectors to m-dimensional ones. In this book matrices are almost always square (n×n) and complex: a gate acting on n qubits is a 2ⁿ×2ⁿ matrix. Element A[i, j] sits at row i, column j — row first, always, in math and in numpy. The matrix is not the operation, it is the operation's representation in a chosen basis; change basis and the same operation gets new entries (4.14). This mirrors software exactly: an interface and its implementation in a given coordinate system. Keep in mind the scale cliff from chapter 3: a two-qubit gate is 4×4 (trivial), a ten-qubit gate is 1024×1024 (16 MB), a thirty-qubit gate is unrepresentable. Matrices are the interface; size is the enemy.
4.2Matrix addition
Matrices of identical shape add elementwise: (A + B)[i, j] = A[i, j] + B[i, j], with the familiar laws — commutative, associative, scalar-distributive, and A + 0 = A. In numpy, A + B and 3*A just work, provided shapes match; mismatched shapes raise errors rather than guessing. Alone, addition is nearly trivial in quantum computing: adding two candidate gate matrices rarely means anything physical. It earns its keep in constructions — Hamiltonians are built as sums of simpler terms (H = H₁ + H₂ + …, each a tensor product, chapter 5), and projector decompositions write an observable as a weighted sum of projectors. The engineering takeaway is small but real: because Hamiltonians are sums of terms, you can often apply each term separately (Lie–Trotter splitting) instead of ever forming the full exponential — a trick the simulation chapters lean on hard.
4.3Matrix multiplication
The product AB is defined when A's column count equals B's row count; entry (AB)[i, j] = Σₖ A[i, k]·B[k, j] — row i of A dotted with column j of B. Two properties separate this from ordinary arithmetic. First, it is generally not commutative: AB ≠ BA, and in quantum mechanics non-commutation is physical — measuring one observable then another is not the same in both orders (the commutator [A, B] = AB − BA quantifies it, and quantum error is largely the story of commuting versus non-commuting terms). Second, composition is the meaning: applying gate B then gate A to a state is the matrix product A·B, right-to-left like function composition — A @ (B @ psi), or precomputed A @ B. The cost is cubic in dimension, O(n³) for n×n — the constant behind "simulation gets slow".
import numpy as np
X = np.array([[0, 1], [1, 0]], dtype=complex)
Z = np.array([[1, 0], [0, -1]], dtype=complex)
print(X @ Z) # anticommuting: XZ = -ZX, the Pauli algebra4.4Matrix-vector multiplication
For square A and vector v, Av is the vector whose i-th component is Σⱼ A[i, j]·v[j] — each output component mixes all input components with weights from row i. This is the fundamental motion of quantum computing: a gate application is Av, updating all 2ⁿ amplitudes at once. In numpy, A @ v (note: A * v broadcasts elementwise and is a common, silent bug). Complexity: O(n²) multiply-adds for dense n×n, O(n²) memory for the matrix itself — for n = 2⁴⁰ qubit-dimension that is already impossible, which is why structured representations (sparse, tensor networks, stabilizers) exist. Two identities worth memorizing: (AB)v = A(Bv) (you can pre-compose gates), and ⟨u|A|v⟩ = ⟨u|(Av) — matrix-vector products and inner products interleave freely, which is how expectation values get computed.
4.5Identity matrices
The identity matrix I has ones on the diagonal and zeros elsewhere: Iv = v for every v. It is the do-nothing operation, the 0 of multiplication, and it anchors every definition in this chapter: an inverse satisfies A·A⁻¹ = I; a unitary matrix satisfies U†U = I; a projector satisfies P² = P. In numpy, np.eye(n). Identities appear in quantum computing in slightly surprising roles: adding a phase with a global factor e^{iφ}·I is physically unobservable (global phase); the decomposition Σᵢ |bᵢ⟩⟨bᵢ| = I (resolution of the identity in an orthonormal basis) is the workhorse step in derivations, letting you insert "a sum over all outcomes" anywhere; and identity on some qubits while others are acted on is the tensor-product construction I ⊗ U (5.11). Do-nothing, done formally, turns out to be a load-bearing object.
4.6Inverse matrices
The inverse A⁻¹ is the matrix with A·A⁻¹ = A⁻¹·A = I; it exists exactly when A is square and full-rank (non-singular). numpy computes it with np.linalg.inv. But the engineering rule is blunt: almost never invert a matrix explicitly. Solving Ax = b as np.linalg.solve(A, b) (LU-based) is faster and numerically safer than computing A⁻¹ then multiplying — explicit inversion can cost an order of magnitude in accuracy for ill-conditioned systems. The deeper reason quantum computing sidesteps inverses: physical evolution is unitary, and unitary matrices are always invertible with U⁻¹ = U† — the conjugate transpose, computable without any factorization. So the general-matrix inverse machinery matters mostly in the classical numerics around your quantum code: fitting calibration models, solving linear systems in state tomography, regularizing least-squares. When you do need one, ask whether a solve would do instead.
4.7Transpose
The transpose Aᵀ flips a matrix across its diagonal: (Aᵀ)[i, j] = A[j, i], turning rows into columns. Algebraic rules: (AB)ᵀ = BᵀAᵀ — order reverses, a pattern that repeats for every transpose-like operation and is worth drilling until automatic. In numpy, A.T. On its own, transpose is a minor player in quantum computing because the amplitudes are complex and conjugation-free transposition misses the essential operation; it appears mainly inside the definition of the conjugate transpose (4.8) and in real-valued classical computations where it costs nothing. One place to meet it early: symmetric real matrices (Aᵀ = A) are the classical cousins of Hermitian matrices, and every theorem you will later use about Hermitian operators — real eigenvalues, orthogonal eigenvectors — has this real-symmetric ancestor. Learn the complex version directly; the real one is the special case.
4.8Conjugate transpose
The conjugate transpose (adjoint, dagger) A† = (Aᵀ)* transposes and conjugates every entry: (A†)[i, j] = A[j, i]*. In numpy, A.conj().T (order does not matter). This is quantum computing's central matrix operation. It defines the adjoint of a vector (the bra), the dagger of a gate (running the gate backwards — uncomputation, the inverse operation every algorithm needs), and the two classes that structure the whole field: unitary matrices (U†U = I, evolution, chapter 4's heart) and Hermitian matrices (H† = H, observables and Hamiltonians). The reversal rule carries over: (AB)† = B†A† — to undo a sequence of gates, undo them in reverse order, exactly like unwinding a call stack. In code, whenever you need "the reverse operation", the answer is almost always .conj().T rather than a numerical inverse.
4.9Hermitian matrices
A matrix H is Hermitian when H† = H — each entry equals the conjugate of its mirror image; diagonal entries must be real. Hermitian matrices are quantum mechanics' real numbers: every observable (energy, spin, parity) is represented by one, because two theorems guarantee physical sanity. First, all eigenvalues are real — measurement outcomes cannot be complex. Second, eigenvectors of distinct eigenvalues are orthogonal — outcomes correspond to distinguishable alternatives. numpy's np.linalg.eigvalsh exploits Hermitian structure, running faster and more stably than the general eig. Two facts the engineer uses weekly: any Hermitian H decomposes as a weighted sum of projectors H = Σ λᵢ Pᵢ (the spectral theorem, 4.18), which is precisely how measurement statistics are computed; and Pauli matrices are the minimal Hermitian alphabet — any 2ⁿ-sized Hermitian is a sum of Pauli tensor products (5.11), the format chemistry Hamiltonians ship in.
4.10Unitary matrices
A matrix U is unitary when U†U = I: it preserves inner products, hence lengths, hence probabilities. Unitarity is quantum mechanics' conservation law — time evolution cannot amplify or discard probability, and because U† = U⁻¹, every quantum operation is freely reversible. The columns (equivalently rows) of a unitary form an orthonormal basis (3.19): unitaries are exactly the transformations mapping one orthonormal basis to another. For the engineer, three consequences dominate. Validation: np.allclose(U.conj().T @ U, np.eye(n)) is the unit test for any gate you construct — run it always. Composition: products and tensor products of unitaries are unitary, so circuits stay valid by construction. Numerics: a simulator that applies only unitaries should keep ‖ψ‖ = 1 forever; if the norm drifts, you have a bug or you have left unitary ground (open systems, noise channels — chapter 33).
4.11Normal matrices
A matrix is normal when A†A = AA† — it commutes with its own adjoint. The class includes all unitary and all Hermitian matrices, plus others (skew-Hermitian, diagonal-with-complex-entries). Its importance is one theorem: a matrix is normal exactly when it is unitarily diagonalizable — there is an orthonormal eigenbasis, A = V·D·V† with V unitary and D diagonal. This is why quantum mechanics can be handled with the clean spectral machinery everywhere: observables (Hermitian) and evolution (unitary) are both normal, both diagonalizable in orthonormal bases, both expressible as "eigen-directions with eigen-values". Non-normal matrices exist in classical numerics (they can be violently transient — growing enormously before decaying) but essentially never as closed-system quantum evolution. The practical check np.allclose(A.conj().T @ A, A @ A.conj().T) tells you whether the cheap, stable eigensolvers apply to your matrix.
4.12Projectors
A projector P satisfies P² = P: applying it twice changes nothing — it projects. The elementary projectors are the rank-1 outer products P = |φ⟩⟨φ| for a unit vector φ, which extract the component of a state along φ and zero out the rest. Hermitian by construction, they are the atoms of quantum measurement: a measurement with outcomes i and associated projectors Pᵢ satisfying Σ Pᵢ = I produces outcome i with probability ⟨ψ|Pᵢ|ψ⟩ and leaves the state Pᵢ|ψ⟩/√(that probability). This is the Born rule in operational form, and it covers computational-basis measurement (P₀ = |0⟩⟨0|, P₁ = |1⟩⟨1|), partial measurement of some qubits (tensor the local projectors with I, 5.12), and observables via the spectral decomposition. In numpy, P = np.outer(phi, phi.conj()). Post-selection — keeping only shots with a chosen outcome — is exactly applying a projector statistically.
4.13Diagonal matrices
A diagonal matrix has nonzeros only on its diagonal: (Dv)[i] = dᵢ·v[i] — each component scaled independently, no mixing. This is the easiest matrix class to compute with: multiplication is O(n) not O(n³), powers are just powers of the entries, functions apply entrywise (4.19), and exponentials are np.exp(d) on the diagonal. Diagonal matrices are also the quantum operations you can actually afford at scale: phase gates, T and S, are diagonal; and any operator is diagonal in its own eigenbasis. The whole strategy of "diagonalize, then act" — change basis to the eigenbasis, scale entrywise, change back — works because conjugation preserves the class structure: V·D·V†. Know also the size trap: a diagonal matrix is still an object of its dimension; at 2ⁿ scale, store it as a length-2ⁿ vector of diagonal entries, never as a full array. Sparse-aware code (7.5) enforces this automatically.
4.14Change of basis
If B is the matrix whose columns are a new orthonormal basis, the coordinates of a vector transform as v_new = B†·v, and operators transform as A_new = B†·A·B — a conjugation. For unitary B this is exact, stable, and reversible: the same vector, same operator, different coordinate description. This single operation explains several quantum-computing staples. The Hadamard gate is the change of basis between the computational basis and the ± basis — apply it, and "amplitudes" become "which-basis-amplitudes". The quantum Fourier transform (chapter 24) is a change of basis into the frequency domain, which is why period-finding sits at the heart of Shor's algorithm. Eigenbasis representation (4.17) is a change of basis that makes operators diagonal. In code: A_new = B.conj().T @ A @ B, with B unitary — verified, as always, by the unitarity check of 4.10.
4.15Eigenvalues
An eigenvalue of A is a scalar λ with a nonzero vector satisfying Av = λv — the transformation stretches that special direction by λ without rotating it. Eigenvalues answer the questions engineers actually ask: what are the possible measurement outcomes (eigenvalues of the observable), at what rates does a system evolve (eigenvalues of the Hamiltonian become phases e^{iλt}), and is this circuit's behavior diagonalizable at all. For Hermitian and unitary matrices — the quantum classes — eigenvalues are real, or unit-magnitude complex respectively: outcomes are real numbers, evolution is pure phase. numpy exposes them via np.linalg.eigvalsh (Hermitian, preferred) and np.linalg.eig (general). The characteristic polynomial det(A − λI) = 0 defines them theoretically, but is numerically useless; production code finds eigenvalues iteratively (7.6). When a chapter later says "the energy levels are the eigenvalues", that sentence is the whole of spectroscopy.
4.16Eigenvectors
Eigenvectors are the directions attached to eigenvalues: the nonzero solutions v of Av = λv, defined only up to scale (and, in quantum computing, fixed by normalization and an arbitrary phase). For a normal matrix, eigenvectors of distinct eigenvalues are orthogonal, so the full set forms an orthonormal eigenbasis — the coordinate system in which the operator becomes diagonal. Physically, eigenvectors are the states that are invariant under the operation: energy eigenstates accumulate only phase under time evolution, and measurement eigenstates are the definite outcomes. Numerically, np.linalg.eigh returns them orthonormalized for Hermitian input. Two habits: degenerate eigenvalues (multiplicity > 1) leave a subspace of valid eigenvectors, any orthonormal basis of which is equally correct — do not compare eigenvectors between runs, compare subspaces; and always check A @ v ≈ λ·v after solving. Eigen-decompositions you have not verified are guesses with formatting.
4.17Eigendecomposition
Eigendecomposition writes a normal matrix as A = V·D·V†: V's columns are the orthonormal eigenvectors, D is diagonal with the eigenvalues. This is the master representation of the chapter because it makes every hard operation easy. Powers: Aᵏ = V·Dᵏ·V†. Exponentials: e^{A} = V·e^{D}·V†, the entrywise exponential of a diagonal matrix — which is exactly how quantum time evolution U(t) = e^{−iHt} is computed for a Hermitian Hamiltonian H, since eigenvalues of H become phases e^{−iλt}. Measurement: the spectral decomposition (4.18) of any observable. In numpy: w, V = np.linalg.eigh(H) (Hermitian) gives you both pieces; reconstruction V @ np.diag(w) @ V.conj().T should reproduce H to roundoff. Cost is O(n³) time and O(n²) memory — fine to a few thousand dimensions, impossible at 2ⁿ for n beyond ~15, which is precisely where chapter 7's iterative methods take over.
import numpy as np
H = np.array([[1, 1j], [-1j, 1]], dtype=complex) # Hermitian
w, V = np.linalg.eigh(H) # eigenvalues, eigenvectors
t = 0.7
U = V @ np.diag(np.exp(-1j * w * t)) @ V.conj().T # U(t) = exp(-i H t)
print(np.allclose(U.conj().T @ U, np.eye(2))) # True: evolution is unitary4.18Spectral decomposition
For any normal A, the spectral decomposition is A = Σᵢ λᵢ·Pᵢ: a weighted sum of projectors Pᵢ onto the eigenspaces — orthogonal, summing to identity. This is the bridge between chapter 3's linear algebra and quantum measurement, and it deserves to be read as an interface specification. Observable H = Σ λᵢ·Pᵢ means: possible outcomes are the λᵢ; probability of outcome i is ⟨ψ|Pᵢ|ψ⟩; post-measurement state is Pᵢ|ψ⟩ normalized. One formula, all of measurement. Functions likewise act spectrally: f(A) = Σ f(λᵢ)·Pᵢ for any scalar function f, which is how you apply a function to an operator without forming matrices (4.19). Degenerate eigenvalues merge into higher-rank projectors — the formula survives unchanged. Numerically, build Pᵢ from eigh's output by grouping eigenvectors with equal eigenvalues; then ⟨ψ|Pᵢ|ψ⟩ is just a couple of inner products — measurement statistics at O(n²) instead of O(n³).
4.19Matrix functions
A matrix function f(A) extends a scalar function to matrices. For diagonalizable A = V·D·V†, the definition is f(A) = V·f(D)·V† — apply f to each eigenvalue. This is not a curiosity: the objects of quantum dynamics are all matrix functions. The time-evolution operator U(t) = e^{−iHt} is the exponential function of the Hamiltonian; e^{A} = Σ Aᵏ/k! converges for every matrix but is numerically treacherous as a naive series — production code uses Padé approximants with scaling-and-squaring (scipy.linalg.expm), or, for Hermitian H with known spectrum, the eigendecomposition route of 4.17. Square roots √A (of positive-definite matrices) appear in fidelity computations; logarithms appear in entropy and in decomposing unitaries into Hamiltonians (U = e^{iH}). The engineering rule: know which route your library takes, because ill-conditioned eigenvalues make one route accurate and another garbage.
4.20Exponentials of matrices
The matrix exponential deserves its own section because e^{A} is the single most consequential matrix function in this book: Schrödinger evolution U(t) = e^{−iHt}, thermal states e^{−H/kT}, and continuous-time quantum walks are all matrix exponentials. Properties: e^{A} is always invertible with inverse e^{−A}; if H is Hermitian, e^{−iHt} is unitary — conservation of probability falls out of the algebra. The product formula e^{A+B} = e^{A}e^{B} holds only when A and B commute; when they do not, the discrepancy — captured by the Baker–Campbell–Hausdorff commutator series — is exactly the Trotter error that quantum simulation algorithms (chapter 30) must bound. Numerically: scipy.linalg.expm for small dense matrices; eigendecomposition when H is Hermitian and its spectrum is cheap; Krylov methods (7.6) for large sparse H. Every quantum simulation, on hardware or on your laptop, is an exercise in computing this object efficiently.