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

Manipulate complex-valued NumPy arrays representing kets and operators.
Use @, np.kron, np.linalg.eig(h), and scipy.linalg.expm fluently for QM.
Verify the algebraic properties we'll demand of operators (unitarity, Hermiticity, normalization).
Plot a single-qubit state on the Bloch sphere and visualize measurement distributions.

All snippets assume:
import numpy as np

1. Complex Numbers and Complex Arrays

Quantum amplitudes are complex. Python's imaginary unit is 1j.

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.

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$
Norm	$\sqrt{\langle\psi	\psi\rangle}$
Outer product	$	\psi\rangle\langle\phi
Matrix–vector	$A	\psi\rangle$
Matrix–matrix	$AB$	`A @ B`

3. Operators as Matrices

The single-qubit gates we'll use constantly:

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

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.

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. 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).

# 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:

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

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

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; here is the computation:

def bloch_vector(psi):
    """Return (x, y, z) = (<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

Why must you use np.vdot (not np.dot) for the inner product of complex kets?
What's the difference between A.T and dag(A), and which do quantum operators need?
Why is np.kron(A, B) order-sensitive, and how does it relate to qubit labeling?
Why is scipy.linalg.expm(M) required for time evolution instead of np.exp(M)?
Two state vectors differ only by a factor 1j. Are they the same physical state? How do you test that numerically?

Answers

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.
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$ ).
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.
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.
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).