# Appendix B — Python/NumPy Refresher for Quantum Computing > You already know Python. This appendix is a **targeted refresher of the NumPy idioms this program > leans on**: complex arrays, the linear-algebra routines we use for states and operators, the > Kronecker product for tensoring qubits, and a couple of plotting helpers (Bloch sphere, state > bar charts) you'll reuse throughout. Skim it, run the snippets, and bookmark it. ## Learning Objectives 1. Manipulate **complex-valued** NumPy arrays representing kets and operators. 2. Use `@`, `np.kron`, `np.linalg.eig(h)`, and `scipy.linalg.expm` fluently for QM. 3. Verify the algebraic properties we'll demand of operators (unitarity, Hermiticity, normalization). 4. Plot a single-qubit state on the **Bloch sphere** and visualize measurement distributions. > All snippets assume: > ```python > import numpy as np > ``` --- ## 1. Complex Numbers and Complex Arrays Quantum amplitudes are complex. Python's imaginary unit is `1j`. ```python z = 2 + 3j z.real, z.imag # (2.0, 3.0) z.conjugate() # (2-3j) abs(z) # 3.605... = sqrt(13) = modulus |z| np.angle(z) # 0.982... = arg(z) in radians # A ket is a complex column vector. dtype=complex is essential. ket0 = np.array([1, 0], dtype=complex) # |0> ket1 = np.array([0, 1], dtype=complex) # |1> plus = (ket0 + ket1) / np.sqrt(2) # |+> = (|0>+|1>)/sqrt(2) ``` > ⚠️ **Pitfall:** `np.array([1, 0])` is integer dtype. Multiplying by `1j` or `np.sqrt` will silently > truncate or warn. **Always** pass `dtype=complex` for states and operators. --- ## 2. Inner Products, Norms, Outer Products For complex vectors the inner product conjugates the left factor: $\langle\psi|\phi\rangle = \sum_i \psi_i^{*}\phi_i$. In NumPy use `np.vdot` (it conjugates its **first** argument), not `np.dot`. ```python np.vdot(plus, plus) # ⟨+|+⟩ = (1+0j) → normalized np.vdot(ket0, ket1) # ⟨0|1⟩ = 0j → orthogonal # Norm of a ket and normalization: def normalize(psi): return psi / np.linalg.norm(psi) # Outer product |ψ⟩⟨ψ| (a projector when |ψ⟩ is normalized): P_plus = np.outer(plus, plus.conj()) # note the explicit conjugate on the bra ``` | Operation | Math | NumPy | |---|---|---| | Inner product | $\langle\psi|\phi\rangle$ | `np.vdot(psi, phi)` | | Norm | $\sqrt{\langle\psi|\psi\rangle}$ | `np.linalg.norm(psi)` | | Outer product | $|\psi\rangle\langle\phi|$ | `np.outer(psi, phi.conj())` | | Matrix–vector | $A|\psi\rangle$ | `A @ psi` | | Matrix–matrix | $AB$ | `A @ B` | --- ## 3. Operators as Matrices The single-qubit gates we'll use constantly: ```python I = np.array([[1, 0], [0, 1]], dtype=complex) X = np.array([[0, 1], [1, 0]], dtype=complex) # Pauli-X (NOT) Y = np.array([[0, -1j], [1j, 0]], dtype=complex) # Pauli-Y Z = np.array([[1, 0], [0, -1]], dtype=complex) # Pauli-Z H = np.array([[1, 1], [1, -1]], dtype=complex) / np.sqrt(2) # Hadamard S = np.array([[1, 0], [0, 1j]], dtype=complex) # phase gate # Adjoint (conjugate transpose) = "dagger": def dag(A): return A.conj().T # Apply a gate to a state: H @ ket0 # = |+> ``` ### Checking the properties we will demand ```python def is_unitary(U, tol=1e-12): return np.allclose(dag(U) @ U, np.eye(U.shape[0]), atol=tol) def is_hermitian(A, tol=1e-12): return np.allclose(A, dag(A), atol=tol) is_unitary(H), is_hermitian(Z), is_hermitian(Y) # (True, True, True) ``` > ⚠️ **Pitfall:** `A.T` is only the transpose. For quantum operators you almost always want the > **conjugate** transpose `A.conj().T`. Define `dag` once and use it everywhere. --- ## 4. Eigenvalues, Spectral Decomposition, Matrix Functions Hermitian operators (observables, Hamiltonians) have real eigenvalues and orthonormal eigenvectors. Use `np.linalg.eigh` for Hermitian matrices — it's faster and returns sorted, real eigenvalues; reserve `np.linalg.eig` for general matrices. ```python vals, vecs = np.linalg.eigh(Z) # vals=[-1, 1]; vecs columns are |1>, |0> # vecs[:, i] is the eigenvector for vals[i]. # Time evolution needs the matrix exponential of a Hamiltonian: U = exp(-i H t). from scipy.linalg import expm def evolve(H, t): return expm(-1j * H * t) # unitary propagator U = evolve(Z, np.pi / 2) is_unitary(U) # True ``` > We derive *why* $U=e^{-iHt}$ in [Term 1.6 · Time Evolution](../term-1-quantum-mechanics/c06-dynamics/01-time-evolution.md). > The matrix exponential is **not** the elementwise `np.exp(M)`; use `scipy.linalg.expm`. --- ## 5. Tensor Products (Combining Qubits) The state of multiple qubits lives in a tensor-product space; the Kronecker product `np.kron` implements it. Ordering matters and must match your convention (we use **big-endian**: qubit 0 is the leftmost/most-significant, matching Braket's bit-string output). ```python # Two-qubit basis state |01> = |0> ⊗ |1> ket01 = np.kron(ket0, ket1) # shape (4,) → [0,1,0,0] # Two-qubit operator: X on qubit 0, I on qubit 1 (X ⊗ I) XI = np.kron(X, I) # shape (4,4) # CNOT (control=0, target=1) built explicitly: CNOT = np.array([[1,0,0,0], [0,1,0,0], [0,0,0,1], [0,0,1,0]], dtype=complex) # Bell state |Φ+> = CNOT (H⊗I) |00> ket00 = np.kron(ket0, ket0) bell = CNOT @ np.kron(H, I) @ ket00 np.round(bell, 3) # [0.707, 0, 0, 0.707] = (|00>+|11>)/sqrt(2) ``` > ⚠️ **Pitfall:** `np.kron(A, B)` ≠ `np.kron(B, A)`. Pick big-endian and stay consistent so your > NumPy results line up with Braket's measurement bit-strings. A small helper for tensoring many factors: ```python from functools import reduce def tensor(*ops): return reduce(np.kron, ops) tensor(H, I, I) # H on qubit 0 of a 3-qubit register ``` --- ## 6. Measurement Probabilities and Sampling Given a state vector, the Born rule gives outcome probabilities $p_i = |\langle i|\psi\rangle|^2 = |\psi_i|^2$. ```python def probabilities(psi): return np.abs(psi) ** 2 def sample(psi, shots=1000, rng=None): rng = rng or np.random.default_rng() p = probabilities(psi) n = int(round(np.log2(len(psi)))) outcomes = rng.choice(len(psi), size=shots, p=p) # Format each integer outcome as an n-bit big-endian string. from collections import Counter return Counter(format(o, f"0{n}b") for o in outcomes) sample(bell, shots=1000) # ≈ Counter({'00': 500, '11': 500}) ``` This mirrors what Braket's `LocalSimulator` does for you; building it once by hand makes the Born rule concrete. --- ## 7. Plotting Helpers ### 7.1 Measurement bar chart ```python import matplotlib.pyplot as plt def plot_counts(counts, title="Measurement outcomes"): keys = sorted(counts) plt.bar(keys, [counts[k] for k in keys]) plt.xlabel("bit-string"); plt.ylabel("count"); plt.title(title) plt.tight_layout(); plt.show() ``` ### 7.2 Bloch-sphere coordinates of a single qubit A pure single-qubit state maps to a point $(\langle X\rangle, \langle Y\rangle, \langle Z\rangle)$ on the unit sphere. We derive this in [Term 1.2 · The Bloch Sphere](../term-1-quantum-mechanics/c02-qubits-bloch/02-bloch-sphere.md); here is the computation: ```python def bloch_vector(psi): """Return (x, y, z) = (, , ) for a single-qubit pure state.""" rho = np.outer(psi, psi.conj()) # density matrix |ψ⟩⟨ψ| x = np.real(np.trace(rho @ X)) y = np.real(np.trace(rho @ Y)) z = np.real(np.trace(rho @ Z)) return np.array([x, y, z]) bloch_vector(ket0) # [0, 0, 1] → north pole bloch_vector(plus) # [1, 0, 0] → +x axis ``` > For a polished 3-D Bloch sphere, the `qutip` library's `Bloch` class is excellent > (`pip install qutip`); we use the lightweight version above to avoid an extra dependency in core > lessons. --- ## 8. Numerical Hygiene - **Tolerances:** floating-point means exact zeros rarely appear. Compare with `np.allclose(a, b, atol=1e-12)`, not `==`. - **Global phase:** $|\psi\rangle$ and $e^{i\theta}|\psi\rangle$ are physically identical. Two state vectors can look different yet be equal up to phase — compare `|⟨φ|ψ⟩|`, not raw components. - **Reproducibility:** seed your RNG with `np.random.default_rng(seed)` so sampled results are repeatable in exercises. - **dtype discipline:** keep states/operators `complex`; mixing in integer or real arrays is the #1 source of silent bugs. --- ## Checkpoint 1. Why must you use `np.vdot` (not `np.dot`) for the inner product of complex kets? 2. What's the difference between `A.T` and `dag(A)`, and which do quantum operators need? 3. Why is `np.kron(A, B)` order-sensitive, and how does it relate to qubit labeling? 4. Why is `scipy.linalg.expm(M)` required for time evolution instead of `np.exp(M)`? 5. Two state vectors differ only by a factor `1j`. Are they the same physical state? How do you test that numerically?
Answers 1. The Hermitian inner product conjugates the bra: $\sum_i \psi_i^* \phi_i$. `np.vdot` conjugates its first argument; `np.dot` does not, giving the wrong value for complex inputs. 2. `A.T` transposes only; `dag(A) = A.conj().T` is the conjugate transpose (adjoint). Quantum operators use the **adjoint** (e.g. unitarity is $U^\dagger U = I$). 3. The Kronecker product is not commutative; `kron(A,B)` places `A` on the more-significant (leftmost, big-endian) qubit and `B` on the next. Swapping them relabels which physical qubit each operator acts on. 4. `expm` computes the true **matrix** exponential $\sum_k M^k/k!$; `np.exp` exponentiates each entry elementwise, which is a different (wrong) object for evolution. 5. Yes — they're identical up to a global phase $e^{i\pi/2}=i$, which is physically unobservable. Test with `np.isclose(abs(np.vdot(a, b)), 1.0)` after normalizing (overlap modulus = 1).
--- ## Further Reading - NumPy linear algebra: - SciPy `expm`: - QuTiP Bloch sphere: - [NC] Nielsen & Chuang, §2.1 (linear algebra) — the math this appendix operationalizes. --- ← Prev: [Appendix A — Braket Setup](A-braket-setup.md) · Next: [Appendix C — Notation & Glossary](C-notation-glossary.md) →