Linear Algebra for Quantum Computing
The essential linear algebra you need for quantum computing: vectors, matrices, tensor products, and inner products, explained with quantum computing in mind.
Quantum computing is built on linear algebra. Every gate, every quantum state, and every measurement has a precise mathematical description in terms of vectors and matrices. This tutorial covers the core ideas you need before studying quantum algorithms or quantum programming frameworks. Every code example uses Python with NumPy, so you can run each block and verify the results yourself.
Vectors as Quantum States
A quantum state is a column vector of complex numbers. For a single qubit, the state lives in a two-dimensional complex vector space. The two basis states, conventionally called |0> and |1>, form the computational basis:
import numpy as np
# The two computational basis states
ket_0 = np.array([1, 0], dtype=complex) # |0>
ket_1 = np.array([0, 1], dtype=complex) # |1>
# A superposition state
alpha = 1 / np.sqrt(2)
beta = 1 / np.sqrt(2)
psi = alpha * ket_0 + beta * ket_1
print(psi) # [0.70710678+0.j 0.70710678+0.j]
The coefficients alpha and beta are called probability amplitudes. They must satisfy the normalization condition: the sum of the squared magnitudes equals 1.
norm_squared = np.sum(np.abs(psi)**2)
print(norm_squared) # 1.0
Column Vectors, Row Vectors, and Dirac Notation
Physicists use Dirac notation (bra-ket notation) as a compact way to write vectors. A ket |psi> is a column vector. A bra <psi| is the conjugate transpose of the corresponding ket, which gives a row vector. In NumPy, taking the conjugate transpose of a 1D array means conjugating each element:
import numpy as np
ket_0 = np.array([1, 0], dtype=complex)
ket_1 = np.array([0, 1], dtype=complex)
# A state with a complex amplitude
psi = np.array([1/np.sqrt(2), 1j/np.sqrt(2)], dtype=complex)
# The bra <psi| is the conjugate transpose
bra_psi = psi.conj()
print("ket |psi> =", psi)
print("bra <psi| =", bra_psi)
# bra <psi| = [0.70710678-0.j 0.-0.70710678j]
The inner product <phi|psi> combines a bra and a ket to produce a single complex number. It measures the overlap between two quantum states:
phi = np.array([1, 0], dtype=complex) # |0>
# Inner product <phi|psi> = phi† @ psi
inner = phi.conj() @ psi
print(f"<phi|psi> = {inner}") # (0.7071067811865476+0j)
The outer product |psi><phi| produces a matrix. When you take the outer product of a state with itself, you get a projector, which projects any vector onto that state:
# Outer product |psi><psi| is a projector (a matrix)
projector = psi[:, None] @ psi[None, :].conj()
print("Projector |psi><psi|:")
print(projector)
# [[0.5+0.j 0.+0.5j]
# [0.-0.5j 0.5+0.j ]]
# Projectors satisfy P^2 = P
print(np.allclose(projector @ projector, projector)) # True
Complex Numbers and Why They Matter
You might wonder whether complex numbers are just a mathematical convenience. They are not. Quantum states require complex amplitudes to describe all physically distinguishable states. Consider two states:
# |+x> = (|0> + |1>) / sqrt(2)
plus_x = np.array([1/np.sqrt(2), 1/np.sqrt(2)], dtype=complex)
# |+y> = (|0> + i|1>) / sqrt(2)
plus_y = np.array([1/np.sqrt(2), 1j/np.sqrt(2)], dtype=complex)
Both states have equal probability of measuring |0> or |1>:
print(f"|+x> P(0)={abs(plus_x[0])**2:.3f}, P(1)={abs(plus_x[1])**2:.3f}")
print(f"|+y> P(0)={abs(plus_y[0])**2:.3f}, P(1)={abs(plus_y[1])**2:.3f}")
# Both give P(0) = 0.500, P(1) = 0.500
Yet these states are physically distinct. They produce different measurement outcomes when measured in other bases (for example, the Y basis). The difference between them lives entirely in the complex phase of the second amplitude. Without complex numbers, quantum computing reduces to classical probabilistic computing; the Gottesman-Knill theorem makes this connection precise.
The Pauli Y gate is one place where complex entries are unavoidable:
Y = np.array([[0, -1j],
[1j, 0]], dtype=complex)
# Y gate applied to |0> gives i|1>
print(Y @ np.array([1, 0], dtype=complex)) # [0.+0.j 0.+1.j]
Matrices as Quantum Gates
A quantum gate is a square matrix that transforms a quantum state. Multiplying a gate matrix by a state vector gives the new state after the gate is applied.
# Hadamard gate: creates superposition from a basis state
H = (1 / np.sqrt(2)) * np.array([[1, 1],
[1, -1]], dtype=complex)
# X gate (NOT gate): flips |0> to |1> and vice versa
X = np.array([[0, 1],
[1, 0]], dtype=complex)
# Apply Hadamard to |0>
result = H @ ket_0
print(result) # [0.70710678+0.j 0.70710678+0.j] -- the |+> state
# Apply X to |0>
result2 = X @ ket_0
print(result2) # [0.+0.j 1.+0.j] -- the |1> state
All valid quantum gates are unitary matrices: a matrix U where U times its conjugate transpose equals the identity. Unitarity guarantees that probabilities sum to 1 after the gate is applied.
def is_unitary(U, tol=1e-10):
return np.allclose(U @ U.conj().T, np.eye(len(U)), atol=tol)
print(is_unitary(H)) # True
print(is_unitary(X)) # True
The Four Pauli Matrices
The Pauli matrices are the most important single-qubit operators. There are four of them, and together they form a basis for all 2x2 Hermitian matrices:
import numpy as np
I = np.eye(2, dtype=complex) # Identity: leaves state unchanged
X = np.array([[0, 1], [1, 0]], dtype=complex) # Bit flip: |0> <-> |1>
Y = np.array([[0, -1j], [1j, 0]], dtype=complex) # Bit + phase flip
Z = np.array([[1, 0], [0, -1]], dtype=complex) # Phase flip: |1> -> -|1>
Each Pauli matrix has a clear physical interpretation:
- I (identity): does nothing. Every state is an eigenvector.
- X (bit flip): swaps |0> and |1>, analogous to a classical NOT gate.
- Y (bit and phase flip): swaps |0> and |1> while also applying complex phases.
- Z (phase flip): leaves |0> unchanged but multiplies |1> by -1.
All four Pauli matrices are both Hermitian (equal to their own conjugate transpose) and unitary (their inverse equals their conjugate transpose). For the non-identity Paulis, this means each matrix is its own inverse:
for name, P in [("I", I), ("X", X), ("Y", Y), ("Z", Z)]:
hermitian = np.allclose(P, P.conj().T)
unitary = np.allclose(P @ P.conj().T, np.eye(2))
print(f"{name}: Hermitian={hermitian}, Unitary={unitary}")
# I: Hermitian=True, Unitary=True
# X: Hermitian=True, Unitary=True
# Y: Hermitian=True, Unitary=True
# Z: Hermitian=True, Unitary=True
Every Pauli matrix has eigenvalues +1 and -1 (the identity has eigenvalue +1 with multiplicity 2):
for name, P in [("X", X), ("Y", Y), ("Z", Z)]:
eigvals = np.linalg.eigvalsh(P)
print(f"{name} eigenvalues: {eigvals}")
# X eigenvalues: [-1. 1.]
# Y eigenvalues: [-1. 1.]
# Z eigenvalues: [-1. 1.]
The fact that {I, X, Y, Z} form a basis for 2x2 Hermitian matrices is why Pauli decompositions appear everywhere in quantum computing. Any observable on a single qubit can be written as a linear combination of these four matrices. This extends to multi-qubit systems: tensor products of Paulis form a basis for all Hermitian matrices on n qubits.
Unitary Matrix Properties
Since every quantum gate must be unitary, it is worth understanding what unitarity buys you. A unitary matrix U satisfies U†U = I, which implies several properties:
import numpy as np
H = (1 / np.sqrt(2)) * np.array([[1, 1], [1, -1]], dtype=complex)
# 1. Columns are orthonormal
col_0, col_1 = H[:, 0], H[:, 1]
print("Column inner products:")
print(f" <c0|c0> = {col_0.conj() @ col_0:.4f}") # 1.0
print(f" <c1|c1> = {col_1.conj() @ col_1:.4f}") # 1.0
print(f" <c0|c1> = {col_0.conj() @ col_1:.4f}") # 0.0
# 2. Rows are orthonormal
row_0, row_1 = H[0, :], H[1, :]
print(f" <r0|r1> = {row_0.conj() @ row_1:.4f}") # 0.0
# 3. Determinant has magnitude 1
det = np.linalg.det(H)
print(f"|det(H)| = {abs(det):.4f}") # 1.0
# 4. Preserves inner products: <phi|psi> = <U phi|U psi>
phi = np.array([1, 0], dtype=complex)
psi = np.array([1/np.sqrt(2), 1j/np.sqrt(2)], dtype=complex)
original_ip = phi.conj() @ psi
transformed_ip = (H @ phi).conj() @ (H @ psi)
print(f"Original inner product: {original_ip:.4f}")
print(f"Transformed inner product: {transformed_ip:.4f}")
print(f"Equal: {np.isclose(original_ip, transformed_ip)}") # True
# 5. Preserves norms: ||U|psi>|| = |||psi>||
norm_before = np.linalg.norm(psi)
norm_after = np.linalg.norm(H @ psi)
print(f"Norm before: {norm_before:.4f}, after: {norm_after:.4f}") # Both 1.0
The preservation of inner products is the mathematical reason that quantum evolution is reversible. Unlike classical logic gates (where AND or OR discard information), unitary gates always preserve the full quantum state, and you can undo any gate by applying its inverse U†.
Inner Products and Probabilities
The inner product of two state vectors gives a complex number that encodes the overlap between states. In quantum mechanics the inner product <phi|psi> is computed by taking the conjugate of phi and dotting it with psi:
import numpy as np
ket_0 = np.array([1, 0], dtype=complex)
ket_1 = np.array([0, 1], dtype=complex)
def inner_product(phi, psi):
return np.dot(phi.conj(), psi)
# Orthogonality of basis states
print(inner_product(ket_0, ket_1)) # 0j
# Normalization
print(inner_product(ket_0, ket_0)) # (1+0j)
The Born rule connects inner products to measurement probabilities. The probability of measuring outcome x from state psi is the squared magnitude of their inner product:
psi = np.array([1/np.sqrt(2), 1/np.sqrt(2)], dtype=complex)
prob_0 = abs(inner_product(ket_0, psi))**2
prob_1 = abs(inner_product(ket_1, psi))**2
print(f"P(0) = {prob_0:.3f}") # 0.500
print(f"P(1) = {prob_1:.3f}") # 0.500
Fidelity Between Quantum States
Fidelity measures how close two quantum states are. For pure states, fidelity is simply the squared magnitude of their inner product:
import numpy as np
ket_0 = np.array([1, 0], dtype=complex)
plus = np.array([1/np.sqrt(2), 1/np.sqrt(2)], dtype=complex)
def fidelity(phi, psi):
"""Fidelity between two pure states: |<phi|psi>|^2"""
return abs(phi.conj() @ psi)**2
# Identical states have fidelity 1
print(f"F(|0>, |0>) = {fidelity(ket_0, ket_0):.3f}") # 1.000
# Orthogonal states have fidelity 0
ket_1 = np.array([0, 1], dtype=complex)
print(f"F(|0>, |1>) = {fidelity(ket_0, ket_1):.3f}") # 0.000
# Partially overlapping states
print(f"F(|0>, |+>) = {fidelity(ket_0, plus):.3f}") # 0.500
# Fidelity between ideal |+> and a noisy approximation
noisy_plus = np.array([0.73, 0.69], dtype=complex)
noisy_plus = noisy_plus / np.linalg.norm(noisy_plus) # normalize
print(f"F(|+>, noisy) = {fidelity(plus, noisy_plus):.4f}")
For density matrices (described later in this tutorial), fidelity generalizes to F(rho, sigma) = (Tr(sqrt(sqrt(rho) sigma sqrt(rho))))^2. The trace distance T(rho, sigma) = (1/2) Tr(|rho - sigma|) provides a complementary measure, where |A| = sqrt(A†A). These two measures are related: 1 - sqrt(F) <= T <= sqrt(1 - F).
Gate Composition: Matrix Multiplication
Applying one gate after another corresponds to matrix multiplication. There is an important subtlety in the ordering: if you apply gate U1 first and then gate U2, the combined operation is U2 @ U1. The last gate applied goes on the left. This is because matrix multiplication acts on the state from the left, so the rightmost matrix hits the state vector first.
import numpy as np
I = np.eye(2, dtype=complex)
H = (1 / np.sqrt(2)) * np.array([[1, 1], [1, -1]], dtype=complex)
X = np.array([[0, 1], [1, 0]], dtype=complex)
ket_0 = np.array([1, 0], dtype=complex)
# Apply H then X: the combined matrix is X @ H
step1 = H @ ket_0 # H|0> = |+>
step2 = X @ step1 # X|+> = |+>
combined = (X @ H) @ ket_0 # Same result in one step
print(np.allclose(step2, combined)) # True
This ordering convention explains why quantum circuit diagrams (which read left to right) look “backwards” compared to matrix expressions (where the rightmost matrix acts first).
Here is a practical example: the CNOT gate can be decomposed as Hadamard on the target, then a controlled-Z, then Hadamard on the target again:
import numpy as np
I = np.eye(2, dtype=complex)
H = (1 / np.sqrt(2)) * np.array([[1, 1], [1, -1]], dtype=complex)
# H on qubit 1 (the target qubit)
H2 = np.kron(I, H)
# Controlled-Z gate
CZ = np.diag([1, 1, 1, -1]).astype(complex)
# CNOT = (I x H) . CZ . (I x H)
CNOT_decomposed = H2 @ CZ @ H2
CNOT_exact = np.array([[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 0, 1],
[0, 0, 1, 0]], dtype=complex)
print(np.allclose(CNOT_decomposed, CNOT_exact)) # True
Tensor Products for Multi-Qubit Systems
When you combine two qubits, the joint state lives in a space that is the tensor product of the two individual spaces. A two-qubit system has dimension 4; an n-qubit system has dimension 2^n. In NumPy, the tensor product (Kronecker product) is computed with np.kron:
import numpy as np
ket_0 = np.array([1, 0], dtype=complex)
ket_1 = np.array([0, 1], dtype=complex)
# Two-qubit basis states via Kronecker product
ket_00 = np.kron(ket_0, ket_0) # [1, 0, 0, 0]
ket_01 = np.kron(ket_0, ket_1) # [0, 1, 0, 0]
ket_10 = np.kron(ket_1, ket_0) # [0, 0, 1, 0]
ket_11 = np.kron(ket_1, ket_1) # [0, 0, 0, 1]
# Product state: qubit 0 in |+>, qubit 1 in |0>
plus = np.array([1/np.sqrt(2), 1/np.sqrt(2)], dtype=complex)
product_state = np.kron(plus, ket_0)
print(product_state) # [0.70710678+0.j 0. +0.j 0.70710678+0.j 0. +0.j]
To apply a single-qubit gate to one qubit in a multi-qubit system, tensor it with the identity on the other qubits:
I = np.eye(2, dtype=complex)
H = (1 / np.sqrt(2)) * np.array([[1, 1], [1, -1]], dtype=complex)
# Hadamard on qubit 0, identity on qubit 1
H_on_q0 = np.kron(H, I)
# Identity on qubit 0, Hadamard on qubit 1
H_on_q1 = np.kron(I, H)
print(H_on_q0.shape) # (4, 4)
The CNOT Gate
The controlled-NOT (CNOT) gate is the most important two-qubit gate. It flips the target qubit if and only if the control qubit is |1>. The full 4x4 matrix is:
import numpy as np
ket_0 = np.array([1, 0], dtype=complex)
ket_1 = np.array([0, 1], dtype=complex)
ket_00 = np.kron(ket_0, ket_0)
ket_01 = np.kron(ket_0, ket_1)
ket_10 = np.kron(ket_1, ket_0)
ket_11 = np.kron(ket_1, ket_1)
CNOT = np.array([[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 0, 1],
[0, 0, 1, 0]], dtype=complex)
# CNOT action on each basis state
for name, state in [('|00>', ket_00), ('|01>', ket_01),
('|10>', ket_10), ('|11>', ket_11)]:
output = CNOT @ state
print(f"CNOT {name} = {output.real}")
# CNOT |00> = [1. 0. 0. 0.] (control=0, no flip)
# CNOT |01> = [0. 1. 0. 0.] (control=0, no flip)
# CNOT |10> = [0. 0. 0. 1.] (control=1, target flipped: |10> -> |11>)
# CNOT |11> = [0. 0. 1. 0.] (control=1, target flipped: |11> -> |10>)
The CNOT gate can create entanglement. Starting from |00> and applying a Hadamard on qubit 0 followed by CNOT produces a Bell state:
I = np.eye(2, dtype=complex)
H = (1 / np.sqrt(2)) * np.array([[1, 1], [1, -1]], dtype=complex)
# Hadamard on qubit 0, identity on qubit 1
H_I = np.kron(H, I)
# Bell state preparation: CNOT @ (H x I) @ |00>
bell = CNOT @ H_I @ ket_00
print("Bell state:", bell)
# [0.70710678+0.j 0.+0.j 0.+0.j 0.70710678+0.j]
# This is (|00> + |11>) / sqrt(2)
Tensor Product Properties
One of the most useful tensor product identities is the mixed product property: (A x B)(C x D) = (AC) x (BD). This means applying A to qubit 0 and B to qubit 1 simultaneously is the same as computing the tensor product of the individual operations:
import numpy as np
H = (1 / np.sqrt(2)) * np.array([[1, 1], [1, -1]], dtype=complex)
X = np.array([[0, 1], [1, 0]], dtype=complex)
Z = np.array([[1, 0], [0, -1]], dtype=complex)
Y = np.array([[0, -1j], [1j, 0]], dtype=complex)
A, B = H, X
C, D = Z, Y
# Mixed product property: (A x B)(C x D) = (AC) x (BD)
left = np.kron(A, B) @ np.kron(C, D)
right = np.kron(A @ C, B @ D)
print(np.allclose(left, right)) # True
This identity is computationally important. If you need to apply single-qubit gates to separate qubits, you can compute each single-qubit product independently (2x2 multiplications) and then take the tensor product, rather than doing full 4x4 or larger matrix multiplications. For n qubits, this can mean the difference between O(n) and O(4^n) operations.
Hermitian Matrices and Observables
An observable is a physical quantity you can measure. In quantum mechanics, observables correspond to Hermitian matrices: matrices equal to their own conjugate transpose. The eigenvalues of a Hermitian matrix are always real, and they are the possible outcomes of measuring that observable.
import numpy as np
Z = np.array([[1, 0],
[0, -1]], dtype=complex)
def is_hermitian(A):
return np.allclose(A, A.conj().T)
print(is_hermitian(Z)) # True
# Eigenvalues of Z are the measurement outcomes: +1 and -1
eigenvalues, eigenvectors = np.linalg.eigh(Z)
print("Outcomes:", eigenvalues) # [-1. 1.]
The expectation value tells you the average measurement outcome for a state psi:
ket_0 = np.array([1, 0], dtype=complex)
ket_1 = np.array([0, 1], dtype=complex)
# Average Z measurement on |+>
psi = np.array([1/np.sqrt(2), 1/np.sqrt(2)], dtype=complex)
expectation = np.real(psi.conj() @ Z @ psi)
print(f"<Z> = {expectation:.3f}") # 0.000 -- equal mix of +1 and -1
Spectral Decomposition
The spectral theorem states that any Hermitian matrix H can be written as H = sum_i lambda_i |e_i><e_i|, where lambda_i are the eigenvalues and |e_i> are the corresponding eigenvectors. This decomposition is central to quantum mechanics: it tells you that measuring an observable H yields outcome lambda_i with probability |<e_i|psi>|^2.
import numpy as np
Z = np.array([[1, 0], [0, -1]], dtype=complex)
eigenvalues, eigenvectors = np.linalg.eigh(Z)
# Reconstruct Z from its spectral decomposition
Z_reconstructed = sum(
lam * np.outer(vec, vec.conj())
for lam, vec in zip(eigenvalues, eigenvectors.T)
)
print(np.allclose(Z, Z_reconstructed)) # True
The spectral decomposition also works for more complex observables. Here is an example with the Pauli X matrix, whose eigenstates are |+> and |->:
X = np.array([[0, 1], [1, 0]], dtype=complex)
eigenvalues, eigenvectors = np.linalg.eigh(X)
print("X eigenvalues:", eigenvalues) # [-1. 1.]
print("X eigenvectors (columns):")
print(eigenvectors)
# The columns are |-> and |+> (up to global phase)
# Measurement probabilities in the X basis
psi = np.array([1, 0], dtype=complex) # |0>
for lam, vec in zip(eigenvalues, eigenvectors.T):
prob = abs(vec.conj() @ psi)**2
print(f"P(X = {lam:+.0f}) = {prob:.3f}")
# P(X = -1) = 0.500
# P(X = +1) = 0.500
Commutators and Pauli Algebra
Two matrices A and B commute if AB = BA. When two gates commute, you can apply them in either order and get the same result. The commutator [A, B] = AB - BA measures how far two operators are from commuting. If [A, B] = 0, the operators commute.
The Pauli matrices have a beautiful algebraic structure. X and Z anti-commute, meaning XZ = -ZX:
import numpy as np
I = np.eye(2, dtype=complex)
X = np.array([[0, 1], [1, 0]], dtype=complex)
Y = np.array([[0, -1j], [1j, 0]], dtype=complex)
Z = np.array([[1, 0], [0, -1]], dtype=complex)
# X and Z anti-commute: XZ = -ZX
print(np.allclose(X @ Z, -(Z @ X))) # True
# The commutator [X, Z] = XZ - ZX = 2 * XZ (since ZX = -XZ)
def commutator(A, B):
return A @ B - B @ A
# Pauli commutation relations: [X, Y] = 2iZ, [Y, Z] = 2iX, [Z, X] = 2iY
print(np.allclose(commutator(X, Y), 2j * Z)) # True
print(np.allclose(commutator(Y, Z), 2j * X)) # True
print(np.allclose(commutator(Z, X), 2j * Y)) # True
These commutation relations matter in several practical contexts:
- Hamiltonian simulation: Trotterization error depends on the commutators of the Hamiltonian terms. If all terms commute, there is no Trotter error at all.
- Quantum error correction: Stabilizer codes are built from commuting sets of Pauli operators. The check for whether an error is detectable involves commutation with these stabilizers.
- Variational algorithms: The gradient of a parameterized quantum circuit involves commutators of the generators with the observable.
The anti-commutator {A, B} = AB + BA is also useful. For distinct Pauli matrices, the anti-commutator is always zero:
def anticommutator(A, B):
return A @ B + B @ A
print(np.allclose(anticommutator(X, Y), np.zeros((2, 2)))) # True
print(np.allclose(anticommutator(Y, Z), np.zeros((2, 2)))) # True
print(np.allclose(anticommutator(Z, X), np.zeros((2, 2)))) # True
Density Matrices
So far, every state we have discussed is a pure state, represented by a single state vector. But quantum systems can also be in mixed states, which arise when there is classical uncertainty about which quantum state a system is in. Density matrices provide a unified framework for both pure and mixed states.
For a pure state |psi>, the density matrix is the outer product rho = |psi><psi|:
import numpy as np
ket_0 = np.array([1, 0], dtype=complex)
ket_1 = np.array([0, 1], dtype=complex)
# Pure state density matrix
psi = np.array([1/np.sqrt(2), 1/np.sqrt(2)], dtype=complex)
rho_pure = np.outer(psi, psi.conj())
print("Pure state density matrix:")
print(rho_pure)
# [[0.5+0.j 0.5+0.j]
# [0.5+0.j 0.5+0.j]]
A mixed state represents classical uncertainty. If a qubit is in state |0> with probability 0.5 and state |1> with probability 0.5, the density matrix is a weighted sum:
# Mixed state: 50% |0>, 50% |1>
rho_mixed = 0.5 * np.outer(ket_0, ket_0.conj()) + 0.5 * np.outer(ket_1, ket_1.conj())
print("Mixed state density matrix:")
print(rho_mixed)
# [[0.5+0.j 0.0+0.j]
# [0.0+0.j 0.5+0.j]]
Notice that the mixed state and the pure state |+> have different density matrices, even though both give 50/50 measurement outcomes in the Z basis. The off-diagonal elements (called coherences) distinguish them. A pure superposition has nonzero coherences; a classical mixture does not.
Every valid density matrix satisfies two properties: trace equal to 1, and positive semidefiniteness. You can check purity with the trace of rho squared:
# Trace is 1 for both
print(f"Tr(rho_pure) = {np.trace(rho_pure).real:.3f}") # 1.000
print(f"Tr(rho_mixed) = {np.trace(rho_mixed).real:.3f}") # 1.000
# Purity: Tr(rho^2) = 1 for pure states, < 1 for mixed
purity_pure = np.trace(rho_pure @ rho_pure).real
purity_mixed = np.trace(rho_mixed @ rho_mixed).real
print(f"Purity of pure state: {purity_pure:.3f}") # 1.000
print(f"Purity of mixed state: {purity_mixed:.3f}") # 0.500
Partial Trace and Entanglement
When two qubits are entangled, neither qubit has a well-defined pure state on its own. The partial trace lets you compute the reduced density matrix of one qubit by “tracing out” the other. If the result is a mixed state, the two qubits are entangled.
import numpy as np
ket_0 = np.array([1, 0], dtype=complex)
ket_1 = np.array([0, 1], dtype=complex)
ket_00 = np.kron(ket_0, ket_0)
ket_11 = np.kron(ket_1, ket_1)
# Bell state: (|00> + |11>) / sqrt(2)
bell = (ket_00 + ket_11) / np.sqrt(2)
rho_bell = np.outer(bell, bell.conj())
# Partial trace over qubit 1 to get the reduced state of qubit 0
# Reshape the 4x4 matrix into a (2,2,2,2) tensor, then trace over qubit 1
rho_q0 = rho_bell.reshape(2, 2, 2, 2).trace(axis1=1, axis2=3)
print("Reduced density matrix of qubit 0:")
print(rho_q0)
# [[0.5+0.j 0.0+0.j]
# [0.0+0.j 0.5+0.j]]
# This is the maximally mixed state -- qubit 0 has no definite state
purity = np.trace(rho_q0 @ rho_q0).real
print(f"Purity: {purity:.3f}") # 0.500
The maximally mixed reduced state confirms that the Bell state is maximally entangled. For a product (unentangled) state, the partial trace would yield a pure state with purity 1.
# Compare with a product state: |+> x |0>
plus = np.array([1/np.sqrt(2), 1/np.sqrt(2)], dtype=complex)
product = np.kron(plus, ket_0)
rho_product = np.outer(product, product.conj())
# Partial trace over qubit 1
rho_q0_product = rho_product.reshape(2, 2, 2, 2).trace(axis1=1, axis2=3)
print("Reduced state of qubit 0 (product state):")
print(rho_q0_product)
# [[0.5+0.j 0.5+0.j]
# [0.5+0.j 0.5+0.j]]
purity_product = np.trace(rho_q0_product @ rho_q0_product).real
print(f"Purity: {purity_product:.3f}") # 1.000 -- pure, as expected
Common Mistakes
Here are pitfalls that frequently trip up people who are new to quantum computing linear algebra:
1. Matrix Multiplication Order
In circuit diagrams, time flows left to right: you draw gate U1 before gate U2. But in the matrix expression, U2 goes on the left: the combined operation is U2 @ U1. This is because the state vector sits on the right, and the first operation to touch it must be the rightmost matrix.
import numpy as np
H = (1 / np.sqrt(2)) * np.array([[1, 1], [1, -1]], dtype=complex)
X = np.array([[0, 1], [1, 0]], dtype=complex)
ket_0 = np.array([1, 0], dtype=complex)
# Circuit: |0> -> H -> X means the matrix expression is X @ H @ |0>
correct = X @ H @ ket_0
wrong = H @ X @ ket_0 # This would be the circuit |0> -> X -> H
print(f"Correct (H then X): {correct}")
print(f"Wrong (X then H): {wrong}")
print(f"Same? {np.allclose(correct, wrong)}") # False
2. Forgetting the Complex Conjugate in Inner Products
The inner product <phi|psi> requires conjugating the first vector. If both vectors have real entries, forgetting the conjugate makes no difference. But with complex amplitudes, it matters:
import numpy as np
phi = np.array([1, 1j], dtype=complex) / np.sqrt(2)
psi = np.array([1, 1], dtype=complex) / np.sqrt(2)
correct_ip = phi.conj() @ psi
wrong_ip = phi @ psi # Missing conjugate!
print(f"Correct <phi|psi> = {correct_ip}") # (0.5-0.5j)
print(f"Wrong phi . psi = {wrong_ip}") # (0.5+0.5j)
3. Using Real dtypes When Complex Numbers Are Needed
NumPy silently discards imaginary parts if you assign a complex value to a real-typed array. Always use dtype=complex for quantum state vectors:
import numpy as np
# Dangerous: real dtype silently drops imaginary parts
bad_state = np.array([1, 0], dtype=float)
Y = np.array([[0, -1j], [1j, 0]], dtype=complex)
result = Y @ bad_state # This works, but...
# If you try to store a complex result in a float array:
container = np.zeros(2, dtype=float)
# container[:] = Y @ np.array([1, 0], dtype=complex) # Would raise error
# Safe: always use complex dtype
good_state = np.array([1, 0], dtype=complex)
result = Y @ good_state
print(result) # [0.+0.j 0.+1.j] -- imaginary part preserved
4. Confusing Tensor Product with Matrix Product
The tensor product (np.kron) and matrix product (@) do very different things. The tensor product combines two systems into a larger one, expanding the dimension. The matrix product applies an operator within the same space, keeping the dimension fixed:
import numpy as np
A = np.array([[1, 0], [0, -1]], dtype=complex) # Z gate
B = np.array([[0, 1], [1, 0]], dtype=complex) # X gate
# Tensor product: 2x2 x 2x2 -> 4x4 (combines two qubit spaces)
tensor = np.kron(A, B)
print(f"Tensor product shape: {tensor.shape}") # (4, 4)
# Matrix product: 2x2 @ 2x2 -> 2x2 (composes operations on one qubit)
product = A @ B
print(f"Matrix product shape: {product.shape}") # (2, 2)
5. Confusing State Norm with Inner Product
The norm of a state vector (which should be 1) is a real number. The inner product of two different states is a complex number. These are related but distinct operations:
import numpy as np
psi = np.array([1/np.sqrt(2), 1j/np.sqrt(2)], dtype=complex)
phi = np.array([1, 0], dtype=complex)
# Norm: always real, always non-negative
norm = np.sqrt((psi.conj() @ psi).real)
print(f"||psi|| = {norm:.4f}") # 1.0000
# Inner product: generally complex
ip = phi.conj() @ psi
print(f"<phi|psi> = {ip}") # (0.7071067811865476+0j)
print(f"|<phi|psi>|^2 = {abs(ip)**2:.4f}") # 0.5000 (measurement probability)
Why This Matters
Every concept in quantum computing maps directly to linear algebra:
- Quantum states are unit vectors in complex vector spaces
- Quantum gates are unitary matrices
- Measurement outcomes are eigenvalues of Hermitian matrices
- Multi-qubit systems are tensor products of individual qubit spaces
- Probabilities come from squared magnitudes of inner products
- Mixed states and decoherence are described by density matrices
- Entanglement is detected through the partial trace
- Gate ordering in circuits maps to matrix multiplication order (reversed)
- Observable decomposition uses the Pauli basis and the spectral theorem
With these foundations in place, you can read the mathematical notation in quantum computing papers and understand what quantum circuit simulators are actually computing under the hood.
Next Steps
- Bra-Ket Notation Explained: the physicist’s shorthand for all of the above
- Complex Numbers in Quantum Computing: why complex amplitudes are essential
- Quantum Gates Explained: a tour of the standard gate set
Was this tutorial helpful?