Quantum Gates as Matrices: The Linear Algebra Behind Quantum Computing
Understand quantum gates as 2x2 unitary matrices, see the explicit matrices for X, Y, Z, H, S, and T gates, and learn how matrix multiplication corresponds to running gates sequentially.
Quantum gates are the quantum equivalent of logic gates in classical computing. But unlike AND, OR, and NOT, which are best understood as truth tables, quantum gates are best understood as matrices. This tutorial builds that understanding from scratch.
Qubits as Vectors
A qubit’s quantum state is a two-dimensional complex vector. The two computational basis states are:
|0⟩ = [1, 0]ᵀ
|1⟩ = [0, 1]ᵀ
A general qubit state is a linear combination (superposition) of these:
|ψ⟩ = α|0⟩ + β|1⟩ = [α, β]ᵀ
where α and β are complex numbers satisfying |α|² + |β|² = 1. This normalization condition ensures that probabilities sum to 1.
import numpy as np
# Computational basis states
ket_0 = np.array([1, 0], dtype=complex)
ket_1 = np.array([0, 1], dtype=complex)
# A superposition state: equal probability of 0 and 1
psi = (1/np.sqrt(2)) * ket_0 + (1/np.sqrt(2)) * ket_1
print(psi) # [0.70710678+0.j 0.70710678+0.j]
print(np.linalg.norm(psi)) # 1.0 (normalized)
Gates as Unitary Matrices
A quantum gate acting on a single qubit is a 2x2 matrix. The key requirement is that the matrix must be unitary: U†U = I, where U† is the conjugate transpose of U.
Unitarity enforces two things:
- It preserves the normalization of quantum states (probabilities still sum to 1 after applying the gate)
- It makes the operation reversible (every quantum gate can be undone)
Applying a gate to a qubit state is just matrix-vector multiplication:
|ψ'⟩ = U|ψ⟩
The Standard Single-Qubit Gates
Pauli X Gate (Quantum NOT)
The X gate flips |0⟩ to |1⟩ and |1⟩ to |0⟩, exactly like a classical NOT gate.
X = [[0, 1],
[1, 0]]
X = np.array([[0, 1],
[1, 0]], dtype=complex)
print(X @ ket_0) # [0.+0.j 1.+0.j] = |1⟩
print(X @ ket_1) # [1.+0.j 0.+0.j] = |0⟩
Pauli Y Gate
Y combines a bit flip and a phase flip:
Y = [[0, -i],
[i, 0]]
Y = np.array([[0, -1j],
[1j, 0]], dtype=complex)
print(Y @ ket_0) # [0.+0.j 0.+1.j]
print(Y @ ket_1) # [0.-1.j 0.+0.j]
Pauli Z Gate (Phase Flip)
Z leaves |0⟩ alone and flips the phase of |1⟩:
Z = [[1, 0],
[0, -1]]
Z = np.array([[1, 0],
[0, -1]], dtype=complex)
print(Z @ ket_0) # [1.+0.j 0.+0.j] = |0⟩ (unchanged)
print(Z @ ket_1) # [ 0.+0.j -1.+0.j] = -|1⟩ (phase flip)
The Pauli gates satisfy the useful relationship X² = Y² = Z² = I.
Hadamard Gate
The Hadamard gate creates equal superpositions. It is arguably the most important single-qubit gate.
H = (1/√2) [[1, 1],
[1, -1]]
H = (1/np.sqrt(2)) * np.array([[1, 1],
[1, -1]], dtype=complex)
plus_state = H @ ket_0
print(plus_state) # [0.70710678+0.j 0.70710678+0.j] = |+⟩
minus_state = H @ ket_1
print(minus_state) # [ 0.70710678+0.j -0.70710678+0.j] = |-⟩
# Hadamard is its own inverse
print(H @ H) # [[1.+0.j 0.+0.j], [0.+0.j 1.+0.j]] = identity
Applying H twice returns the original state, which confirms H is unitary and self-inverse (H† = H).
S Gate (Phase Gate)
The S gate applies a 90-degree (π/2) phase rotation to |1⟩:
S = [[1, 0],
[0, i]]
S = np.array([[1, 0],
[0, 1j]], dtype=complex)
print(S @ ket_1) # [0.+0.j 0.+1.j] (phase multiplied by i)
print(S @ S) # This equals Z
T Gate (π/8 Gate)
The T gate applies a 45-degree (π/4) phase rotation:
T = [[1, 0],
[0, e^(iπ/4)]]
T = np.array([[1, 0],
[0, np.exp(1j * np.pi / 4)]], dtype=complex)
print(T @ T) # This equals S (two T gates = one S gate)
The T gate is important in fault-tolerant quantum computing because it is “magic”; it cannot be implemented exactly on many hardware platforms and must be approximated using a technique called magic state distillation.
Verifying Unitarity
You can verify that any gate U satisfies U†U = I:
def is_unitary(U, tol=1e-10):
product = U.conj().T @ U
return np.allclose(product, np.eye(len(U)), atol=tol)
for name, gate in [("X", X), ("Y", Y), ("Z", Z), ("H", H), ("S", S), ("T", T)]:
print(f"{name}: unitary = {is_unitary(gate)}")
# All print True
Eigenvalues and Eigenstates of the Pauli Gates
Every quantum gate has eigenstates: states that, when the gate acts on them, come back multiplied by a scalar (the eigenvalue). Understanding eigenstates is essential because they are the states that remain unchanged (up to a phase) under the corresponding gate operation. Physically, if you measure in the basis defined by a gate’s eigenstates, those eigenstates are the outcomes that survive the measurement process undisturbed.
X Gate Eigenstates
The X gate has eigenvalues +1 and -1. The corresponding eigenstates are:
|+⟩ = (|0⟩ + |1⟩) / √2 with eigenvalue +1
|-⟩ = (|0⟩ - |1⟩) / √2 with eigenvalue -1
This means X|+⟩ = +|+⟩ and X|-⟩ = -|-⟩. Applying X to |+⟩ leaves it completely unchanged. Applying X to |-⟩ flips its global sign, which is physically unobservable, so |-⟩ is also “stable” under X in a measurement sense.
Z Gate Eigenstates
The Z gate has eigenvalues +1 and -1, with eigenstates that are simply the computational basis states:
|0⟩ with eigenvalue +1
|1⟩ with eigenvalue -1
This makes sense: Z|0⟩ = |0⟩ and Z|1⟩ = -|1⟩.
Y Gate Eigenstates
The Y gate has eigenvalues +1 and -1, with eigenstates:
|+i⟩ = (|0⟩ + i|1⟩) / √2 with eigenvalue +1
|-i⟩ = (|0⟩ - i|1⟩) / √2 with eigenvalue -1
Let’s verify all of this numerically using numpy.linalg.eig:
# Compute eigenvalues and eigenvectors for all three Pauli gates
for name, gate in [("X", X), ("Y", Y), ("Z", Z)]:
eigenvalues, eigenvectors = np.linalg.eig(gate)
print(f"\n{name} gate:")
print(f" Eigenvalues: {eigenvalues}")
for i, val in enumerate(eigenvalues):
vec = eigenvectors[:, i]
# Normalize so the first nonzero component is positive real
vec = vec / vec[np.argmax(np.abs(vec))]
vec = vec / np.linalg.norm(vec)
print(f" Eigenvalue {val:+.1f}: eigenvector {vec}")
# Manually verify X|+⟩ = |+⟩ and X|-⟩ = -|-⟩
plus = (ket_0 + ket_1) / np.sqrt(2)
minus = (ket_0 - ket_1) / np.sqrt(2)
print("\nVerification:")
print(f"X|+⟩ = {X @ plus}") # [0.707, 0.707] = |+⟩
print(f"|+⟩ = {plus}") # [0.707, 0.707]
print(f"X|-⟩ = {X @ minus}") # [-0.707, 0.707] = -|-⟩
print(f"-|-⟩ = {-minus}") # [-0.707, 0.707]
print(f"X|+⟩ == |+⟩: {np.allclose(X @ plus, plus)}") # True
print(f"X|-⟩ == -|-⟩: {np.allclose(X @ minus, -minus)}") # True
# Verify Y eigenstates
plus_i = (ket_0 + 1j * ket_1) / np.sqrt(2)
minus_i = (ket_0 - 1j * ket_1) / np.sqrt(2)
print(f"\nY|+i⟩ == |+i⟩: {np.allclose(Y @ plus_i, plus_i)}") # True
print(f"Y|-i⟩ == -|-i⟩: {np.allclose(Y @ minus_i, -minus_i)}") # True
The eigenstates of X form the X-basis (also called the Hadamard basis), the eigenstates of Z form the Z-basis (the computational basis), and the eigenstates of Y form the Y-basis. These three bases are mutually unbiased: a state that is definite in one basis is maximally uncertain in the other two.
Global Phase vs. Relative Phase
Two quantum states that differ only by a global phase factor are physically identical. If |ψ⟩ is a quantum state, then e^{iφ}|ψ⟩ produces the exact same measurement statistics for every possible measurement. No experiment can distinguish them.
However, relative phase between components of a superposition is physically meaningful. Consider the states |+⟩ = (|0⟩ + |1⟩)/√2 and |-⟩ = (|0⟩ - |1⟩)/√2. In the Z basis, both give 50% probability for |0⟩ and 50% probability for |1⟩. They look identical if you only measure in the Z basis. But in the X basis, they are completely distinguishable: |+⟩ is 100% the +1 eigenstate of X, and |-⟩ is 100% the -1 eigenstate of X.
# Define the states
plus = (ket_0 + ket_1) / np.sqrt(2) # |+⟩
minus = (ket_0 - ket_1) / np.sqrt(2) # |-⟩
# Z-basis measurement probabilities: P(0) = |<0|psi>|^2, P(1) = |<1|psi>|^2
print("Z-basis probabilities:")
print(f" |+⟩: P(0) = {abs(ket_0 @ plus)**2:.3f}, P(1) = {abs(ket_1 @ plus)**2:.3f}")
print(f" |-⟩: P(0) = {abs(ket_0 @ minus)**2:.3f}, P(1) = {abs(ket_1 @ minus)**2:.3f}")
# Both give 0.500, 0.500
# X-basis measurement probabilities: project onto |+⟩ and |-⟩
print("\nX-basis probabilities:")
print(f" |+⟩: P(+) = {abs(plus @ plus)**2:.3f}, P(-) = {abs(minus @ plus)**2:.3f}")
print(f" |-⟩: P(+) = {abs(plus @ minus)**2:.3f}, P(-) = {abs(minus @ minus)**2:.3f}")
# |+⟩ gives P(+)=1.000, P(-)=0.000
# |-⟩ gives P(+)=0.000, P(-)=1.000
# Demonstrate that global phase is unobservable
psi = (ket_0 + 1j * ket_1) / np.sqrt(2)
psi_with_phase = np.exp(1j * 0.7) * psi # global phase of 0.7 radians
print("\nGlobal phase check:")
print(f" P(0) original: {abs(ket_0 @ psi)**2:.6f}")
print(f" P(0) with phase: {abs(ket_0 @ psi_with_phase)**2:.6f}")
print(f" P(1) original: {abs(ket_1 @ psi)**2:.6f}")
print(f" P(1) with phase: {abs(ket_1 @ psi_with_phase)**2:.6f}")
# Identical probabilities in every case
This distinction is critical for understanding interference. Quantum algorithms work by arranging relative phases so that correct answers constructively interfere (amplitudes add) and incorrect answers destructively interfere (amplitudes cancel). If relative phase were unobservable like global phase, quantum computing would not offer any advantage over classical computing.
Rotation Gates
The Pauli gates are special cases of a more general family: the rotation gates. These gates rotate the qubit state around the X, Y, or Z axis of the Bloch sphere by a specified angle θ.
Rotation as Matrix Exponential
Each rotation gate is the matrix exponential of the corresponding Pauli operator:
Rx(θ) = exp(-iθX/2) = I·cos(θ/2) - iX·sin(θ/2)
Ry(θ) = exp(-iθY/2) = I·cos(θ/2) - iY·sin(θ/2)
Rz(θ) = exp(-iθZ/2) = I·cos(θ/2) - iZ·sin(θ/2)
The explicit matrices are:
Rx(θ) = [[cos(θ/2), -i·sin(θ/2)],
[-i·sin(θ/2), cos(θ/2) ]]
Ry(θ) = [[cos(θ/2), -sin(θ/2)],
[sin(θ/2), cos(θ/2)]]
Rz(θ) = [[e^(-iθ/2), 0 ],
[0, e^(iθ/2)]]
I = np.eye(2, dtype=complex)
def Rx(theta):
"""Rotation around X axis by angle theta."""
return np.cos(theta/2) * I - 1j * np.sin(theta/2) * X
def Ry(theta):
"""Rotation around Y axis by angle theta."""
return np.cos(theta/2) * I - 1j * np.sin(theta/2) * Y
def Rz(theta):
"""Rotation around Z axis by angle theta."""
return np.cos(theta/2) * I - 1j * np.sin(theta/2) * Z
# Verify the explicit matrix form of Rx
theta = 0.7
Rx_formula = Rx(theta)
Rx_explicit = np.array([
[np.cos(theta/2), -1j * np.sin(theta/2)],
[-1j * np.sin(theta/2), np.cos(theta/2)]
], dtype=complex)
print(f"Rx formula matches explicit: {np.allclose(Rx_formula, Rx_explicit)}") # True
# Verify the explicit matrix form of Rz
Rz_formula = Rz(theta)
Rz_explicit = np.array([
[np.exp(-1j * theta/2), 0],
[0, np.exp(1j * theta/2)]
], dtype=complex)
print(f"Rz formula matches explicit: {np.allclose(Rz_formula, Rz_explicit)}") # True
Special Cases
At specific angles, the rotation gates reduce to familiar gates (up to global phase):
# Rx(π) = -iX
print(f"Rx(π) == -iX: {np.allclose(Rx(np.pi), -1j * X)}") # True
# Rz(π) = -iZ
print(f"Rz(π) == -iZ: {np.allclose(Rz(np.pi), -1j * Z)}") # True
# Ry(π) = -iY
print(f"Ry(π) == -iY: {np.allclose(Ry(np.pi), -1j * Y)}") # True
# The Hadamard gate is Ry(π/2) followed by Rz(π), up to global phase.
# H = Rz(π) @ Ry(π/2), up to global phase
product = Rz(np.pi) @ Ry(np.pi / 2)
# Find the global phase factor
phase = H[0, 0] / product[0, 0]
print(f"H == phase * Rz(π)Ry(π/2): {np.allclose(H, phase * product)}") # True
print(f"Global phase factor: {phase}")
This tells you that the Hadamard gate is not a “special” gate requiring its own hardware implementation. It can be decomposed into rotations, which are the native operations on most quantum hardware.
Gate Decompositions: Euler Angles
Any single-qubit unitary gate U in SU(2) can be decomposed into three rotations:
U = e^{iα} · Rz(γ) · Ry(β) · Rz(δ)
for some angles α (global phase), β, γ, and δ. This is the Euler angle decomposition (also called the ZYZ decomposition). It means that if your hardware can perform Ry and Rz rotations, you can implement any single-qubit gate.
An equivalent decomposition uses Rz and Rx instead:
U = e^{iα} · Rz(γ) · Rx(β) · Rz(δ)
Here is how to extract the Euler angles from an arbitrary 2x2 unitary and reconstruct it:
def decompose_zyz(U):
"""Decompose a 2x2 unitary into Rz(gamma) @ Ry(beta) @ Rz(delta)
plus a global phase.
Returns (global_phase, gamma, beta, delta).
"""
# Ensure U is unitary
assert np.allclose(U.conj().T @ U, np.eye(2), atol=1e-10), "Not unitary"
# Extract the overall phase so det(U') = 1
det_U = np.linalg.det(U)
global_phase = np.angle(det_U) / 2
U_su2 = U * np.exp(-1j * global_phase)
# For a matrix in SU(2):
# U_su2 = [[a, -b*], [b, a*]] with |a|^2 + |b|^2 = 1
a = U_su2[0, 0]
b = U_su2[1, 0]
# beta from |a| = cos(beta/2)
beta = 2 * np.arccos(np.clip(np.abs(a), 0, 1))
if np.abs(beta) < 1e-10:
# beta ~ 0: only Rz matters, split arbitrarily
gamma = np.angle(a)
delta = np.angle(a)
gamma = gamma + delta
delta = 0.0
elif np.abs(beta - np.pi) < 1e-10:
# beta ~ pi: cos(beta/2) ~ 0
gamma = np.angle(b)
delta = -np.angle(b)
gamma = gamma + delta
delta = 0.0
else:
gamma = np.angle(a) + np.angle(b)
delta = np.angle(a) - np.angle(b)
return global_phase, gamma, beta, delta
def reconstruct_zyz(global_phase, gamma, beta, delta):
"""Reconstruct U from Euler angles."""
return np.exp(1j * global_phase) * Rz(gamma) @ Ry(beta) @ Rz(delta)
# Test with the Hadamard gate
gp, gamma, beta, delta = decompose_zyz(H)
H_reconstructed = reconstruct_zyz(gp, gamma, beta, delta)
print(f"H decomposition angles: phase={gp:.4f}, gamma={gamma:.4f}, "
f"beta={beta:.4f}, delta={delta:.4f}")
print(f"Reconstruction matches H: {np.allclose(H_reconstructed, H)}") # True
# Test with the S gate
gp, gamma, beta, delta = decompose_zyz(S)
S_reconstructed = reconstruct_zyz(gp, gamma, beta, delta)
print(f"S reconstruction matches: {np.allclose(S_reconstructed, S)}") # True
# Test with a random unitary
from scipy.stats import unitary_group
U_random = unitary_group.rvs(2)
gp, gamma, beta, delta = decompose_zyz(U_random)
U_recon = reconstruct_zyz(gp, gamma, beta, delta)
print(f"Random U reconstruction matches: {np.allclose(U_recon, U_random)}") # True
Commutators and Anti-Commutators
Two matrices A and B commute if AB = BA, meaning the order of application does not matter. In quantum mechanics, non-commuting operators represent observables that cannot be simultaneously determined with arbitrary precision (this is the generalized uncertainty principle).
The Pauli gates have elegant commutation relations. The commutator [A, B] = AB - BA measures how much two operators fail to commute:
[X, Y] = XY - YX = 2iZ
[Y, Z] = YZ - ZY = 2iX
[Z, X] = ZX - XZ = 2iY
The anti-commutator {A, B} = AB + BA for the Pauli gates is:
{X, Y} = XY + YX = 0
{Y, Z} = YZ + ZY = 0
{Z, X} = ZX + XZ = 0
And each Pauli gate anti-commutes with itself to give twice the identity:
{X, X} = 2I, {Y, Y} = 2I, {Z, Z} = 2I
# Verify commutation relations
print("Commutators:")
print(f" [X, Y] = 2iZ: {np.allclose(X @ Y - Y @ X, 2j * Z)}") # True
print(f" [Y, Z] = 2iX: {np.allclose(Y @ Z - Z @ Y, 2j * X)}") # True
print(f" [Z, X] = 2iY: {np.allclose(Z @ X - X @ Z, 2j * Y)}") # True
print("\nAnti-commutators:")
print(f" {{X, Y}} = 0: {np.allclose(X @ Y + Y @ X, np.zeros((2,2)))}") # True
print(f" {{Y, Z}} = 0: {np.allclose(Y @ Z + Z @ Y, np.zeros((2,2)))}") # True
print(f" {{Z, X}} = 0: {np.allclose(Z @ X + X @ Z, np.zeros((2,2)))}") # True
# Demonstrate that non-commuting gates change results when reordered
state = (ket_0 + 1j * ket_1) / np.sqrt(2)
result_XY = Y @ X @ state # Apply X then Y
result_YX = X @ Y @ state # Apply Y then X
print(f"\nX then Y: {result_XY}")
print(f"Y then X: {result_YX}")
print(f"Same result? {np.allclose(result_XY, result_YX)}") # False
print(f"Difference: {result_XY - result_YX}")
The practical implication: you cannot freely reorder gates in a quantum circuit. Swapping the order of non-commuting gates changes the output state. When optimizing circuits, you can only rearrange gates that commute with each other.
Measurement as a Matrix Operation
In quantum computing, measurement is not just “looking at the qubit.” It is a mathematically precise operation described by projection operators.
Measuring a single qubit in the Z basis (the computational basis) corresponds to projecting the state onto either |0⟩ or |1⟩. The two projection operators are:
P₀ = |0⟩⟨0| = [[1, 0], P₁ = |1⟩⟨1| = [[0, 0],
[0, 0]] [0, 1]]
These projectors satisfy three important properties: they are Hermitian (P† = P), they are idempotent (P² = P), and they form a complete set (P₀ + P₁ = I).
For a state |ψ⟩ = α|0⟩ + β|1⟩:
- Probability of outcome 0: P(0) = ⟨ψ|P₀|ψ⟩ = |α|²
- Probability of outcome 1: P(1) = ⟨ψ|P₁|ψ⟩ = |β|²
- Post-measurement state after getting 0: P₀|ψ⟩ / √P(0) = |0⟩
- Post-measurement state after getting 1: P₁|ψ⟩ / √P(1) = |1⟩
# Define projectors
P0 = np.array([[1, 0], [0, 0]], dtype=complex) # |0><0|
P1 = np.array([[0, 0], [0, 1]], dtype=complex) # |1><1|
def projective_measurement(state, outcome):
"""Compute the probability and post-measurement state for a Z-basis measurement.
Args:
state: a 2-element complex numpy array (qubit state vector)
outcome: 0 or 1
Returns:
(probability, post_measurement_state)
"""
P = P0 if outcome == 0 else P1
projected = P @ state
probability = np.real(state.conj() @ P @ state)
if probability < 1e-15:
return probability, None # outcome is impossible
post_state = projected / np.sqrt(probability)
return probability, post_state
# Test with |+⟩ = (|0⟩ + |1⟩) / √2
plus = (ket_0 + ket_1) / np.sqrt(2)
prob_0, post_0 = projective_measurement(plus, 0)
prob_1, post_1 = projective_measurement(plus, 1)
print(f"|+⟩ measurement:")
print(f" P(0) = {prob_0:.3f}, post-state = {post_0}") # 0.500, [1, 0]
print(f" P(1) = {prob_1:.3f}, post-state = {post_1}") # 0.500, [0, 1]
# Test with a biased state
psi = np.array([np.sqrt(0.3), np.sqrt(0.7)], dtype=complex)
prob_0, post_0 = projective_measurement(psi, 0)
prob_1, post_1 = projective_measurement(psi, 1)
print(f"\nBiased state measurement:")
print(f" P(0) = {prob_0:.3f}, post-state = {post_0}") # 0.300, [1, 0]
print(f" P(1) = {prob_1:.3f}, post-state = {post_1}") # 0.700, [0, 1]
print(f" P(0) + P(1) = {prob_0 + prob_1:.3f}") # 1.000
Measurement collapses the superposition. Before measurement, the state can have nonzero amplitudes for both |0⟩ and |1⟩. After measurement, it is in a definite basis state. This is irreversible, which is why measurement is the one non-unitary operation in quantum computing.
Sequential Gates as Matrix Multiplication
Applying gate A then gate B is the same as multiplying the matrices B @ A. Note the order reversal: the last gate applied goes on the left.
# Apply H then X
# This is X applied after H: result = X @ H @ state
state = ket_0
after_H = H @ state
after_X_H = X @ after_H
# Equivalent: compose gates first
composed = X @ H
print(np.allclose(composed @ ket_0, after_X_H)) # True
Two-Qubit States and Tensor Products
A two-qubit system lives in a 4-dimensional space. The basis states are |00⟩, |01⟩, |10⟩, |11⟩. A product state like |0⟩⊗|1⟩ (qubit 0 in state |0⟩ and qubit 1 in state |1⟩) is computed via the tensor product:
# |01⟩ = |0⟩ ⊗ |1⟩
state_01 = np.kron(ket_0, ket_1)
print(state_01) # [0.+0.j 1.+0.j 0.+0.j 0.+0.j]
# |10⟩ = |1⟩ ⊗ |0⟩
state_10 = np.kron(ket_1, ket_0)
print(state_10) # [0.+0.j 0.+0.j 1.+0.j 0.+0.j]
To apply a gate to only one qubit in a two-qubit system, take the tensor product of that gate with the identity:
# Apply H to qubit 0, leave qubit 1 alone
H_on_qubit0 = np.kron(H, I) # 4x4 matrix
print(H_on_qubit0 @ state_01)
The CNOT Gate as a 4x4 Matrix
The CNOT (Controlled-NOT) gate flips the target qubit if and only if the control qubit is |1⟩. Its matrix is:
CNOT = [[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 0, 1],
[0, 0, 1, 0]]
CNOT = np.array([[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 0, 1],
[0, 0, 1, 0]], dtype=complex)
# |10⟩ -> |11⟩ (control is 1, so target flips from 0 to 1)
print(CNOT @ state_10) # [0.+0.j 0.+0.j 0.+0.j 1.+0.j] = |11⟩
# |00⟩ -> |00⟩ (control is 0, no flip)
state_00 = np.kron(ket_0, ket_0)
print(CNOT @ state_00) # [1.+0.j 0.+0.j 0.+0.j 0.+0.j] = |00⟩
Controlled Gates in Full Generality
The CNOT gate is a specific example of a controlled gate: it applies the X gate to the target qubit when the control qubit is |1⟩. You can build a controlled version of any single-qubit gate G. The general formula is:
C(G) = |0⟩⟨0| ⊗ I + |1⟩⟨1| ⊗ G
In matrix form, this is a 4x4 block matrix:
C(G) = [[1, 0, 0, 0 ],
[0, 1, 0, 0 ],
[0, 0, G[0,0], G[0,1]],
[0, 0, G[1,0], G[1,1]]]
When the control qubit is |0⟩, nothing happens. When it is |1⟩, the gate G acts on the target.
def controlled_gate(G):
"""Build the controlled version of a single-qubit gate G.
The control is qubit 0 (leftmost/most significant),
the target is qubit 1.
"""
proj_0 = np.outer(ket_0, ket_0) # |0><0|
proj_1 = np.outer(ket_1, ket_1) # |1><1|
return np.kron(proj_0, I) + np.kron(proj_1, G)
# Build CZ (controlled-Z)
CZ = controlled_gate(Z)
print("CZ gate:")
print(np.round(CZ, 3))
# [[1, 0, 0, 0],
# [0, 1, 0, 0],
# [0, 0, 1, 0],
# [0, 0, 0, -1]]
# Build controlled-S
CS = controlled_gate(S)
print("\nControlled-S gate:")
print(np.round(CS, 3))
# Build controlled-T
CT = controlled_gate(T)
print("\nControlled-T gate:")
print(np.round(CT, 3))
# Verify the well-known identity: CZ = (H ⊗ I) @ CNOT @ (H ⊗ I)
# This says you can build CZ from CNOT by conjugating the target with H
HI = np.kron(H, I)
# Note: H is applied to qubit 0, but for CZ = (I⊗H)CNOT(I⊗H) on the target
IH = np.kron(I, H)
CZ_from_CNOT = IH @ CNOT @ IH
print(f"\nCZ == (I⊗H) CNOT (I⊗H): {np.allclose(CZ, CZ_from_CNOT)}") # True
# Test CZ on |11⟩: should give -|11⟩
state_11 = np.kron(ket_1, ket_1)
result = CZ @ state_11
print(f"CZ|11⟩ = {result}") # [0, 0, 0, -1] = -|11⟩
The Toffoli (CCX) Gate
The Toffoli gate is a three-qubit gate that flips the target qubit if and only if both control qubits are |1⟩. It is the quantum version of a classical AND gate (in a reversible form). Its 8x8 matrix is the identity everywhere except the bottom-right 2x2 block, which contains the X gate:
Toffoli = diag(1, 1, 1, 1, 1, 1) ⊕ X
In the computational basis ordered as |000⟩, |001⟩, …, |111⟩, the Toffoli gate maps |110⟩ to |111⟩ and |111⟩ to |110⟩, leaving all other basis states unchanged.
# Build the 8x8 Toffoli gate
Toffoli = np.eye(8, dtype=complex)
# Swap the last two rows/columns (the |110⟩ and |111⟩ subspace)
Toffoli[6, 6] = 0
Toffoli[6, 7] = 1
Toffoli[7, 6] = 1
Toffoli[7, 7] = 0
print("Toffoli gate:")
print(Toffoli.real.astype(int))
# Basis states for three qubits
def ket3(a, b, c):
"""Create a 3-qubit computational basis state |abc⟩."""
return np.kron(np.kron(
ket_0 if a == 0 else ket_1,
ket_0 if b == 0 else ket_1),
ket_0 if c == 0 else ket_1)
# Test: |110⟩ -> |111⟩ (both controls are 1, so target flips)
result = Toffoli @ ket3(1, 1, 0)
print(f"\nToffoli|110⟩ = |111⟩: {np.allclose(result, ket3(1, 1, 1))}") # True
# Test: |010⟩ -> |010⟩ (first control is 0, so nothing happens)
result = Toffoli @ ket3(0, 1, 0)
print(f"Toffoli|010⟩ = |010⟩: {np.allclose(result, ket3(0, 1, 0))}") # True
# Test: |100⟩ -> |100⟩ (second control is 0, so nothing happens)
result = Toffoli @ ket3(1, 0, 0)
print(f"Toffoli|100⟩ = |100⟩: {np.allclose(result, ket3(1, 0, 0))}") # True
# Verify unitarity
print(f"Toffoli is unitary: {is_unitary(Toffoli)}") # True
The Toffoli gate is classically universal: you can build any classical Boolean function using Toffoli gates alone. Combined with the Hadamard gate, the Toffoli gate becomes quantum universal, meaning the set {Toffoli, H} can approximate any unitary transformation on any number of qubits to arbitrary precision.
The SWAP Gate and Multi-Qubit Composition
The SWAP gate exchanges the states of two qubits: SWAP|ab⟩ = |ba⟩. Its matrix in the |00⟩, |01⟩, |10⟩, |11⟩ basis is:
SWAP = [[1, 0, 0, 0],
[0, 0, 1, 0],
[0, 1, 0, 0],
[0, 0, 0, 1]]
A key result in quantum computing is that the SWAP gate can be decomposed into exactly three CNOT gates:
SWAP = CNOT(a,b) @ CNOT(b,a) @ CNOT(a,b)
where CNOT(a,b) has qubit a as control and qubit b as target.
SWAP = np.array([[1, 0, 0, 0],
[0, 0, 1, 0],
[0, 1, 0, 0],
[0, 0, 0, 1]], dtype=complex)
# Verify SWAP|10⟩ = |01⟩
print(f"SWAP|10⟩ = |01⟩: {np.allclose(SWAP @ state_10, state_01)}") # True
print(f"SWAP|01⟩ = |10⟩: {np.allclose(SWAP @ state_01, state_10)}") # True
# Build CNOT with reversed control/target (qubit 1 controls, qubit 0 is target)
# CNOT_reverse = (H ⊗ H) @ CNOT @ (H ⊗ H)
# Or build it directly: |0><0| ⊗ I + |1><1| ⊗ X for standard,
# and I ⊗ |0><0| + X ⊗ |1><1| for reversed
CNOT_10 = CNOT # control=qubit 0, target=qubit 1
CNOT_01 = np.kron(I, np.outer(ket_0, ket_0)) + np.kron(X, np.outer(ket_1, ket_1))
# CNOT_01 has control=qubit 1, target=qubit 0
# Verify SWAP = three CNOTs
SWAP_from_CNOTs = CNOT_10 @ CNOT_01 @ CNOT_10
print(f"SWAP from 3 CNOTs: {np.allclose(SWAP, SWAP_from_CNOTs)}") # True
# Build CSWAP (Fredkin gate) using Toffoli and CNOTs
# CSWAP = controlled-SWAP: swap qubits 1 and 2 if qubit 0 is |1⟩
# One construction: Toffoli(0,2,1) @ Toffoli(0,1,2) @ Toffoli(0,2,1)
# But a simpler one uses CNOT + Toffoli:
# Step 1: CNOT with control=qubit2, target=qubit1
# Step 2: Toffoli with controls=qubit0,qubit1, target=qubit2
# Step 3: CNOT with control=qubit2, target=qubit1
# For simplicity, build CSWAP directly and verify
CSWAP = np.eye(8, dtype=complex)
# CSWAP swaps qubits 1 and 2 when qubit 0 is |1⟩
# |100⟩ <-> |100⟩ (no swap needed, both targets are 0)
# |101⟩ <-> |110⟩ (swap the last two qubits)
# The affected basis states are |101⟩ (index 5) and |110⟩ (index 6)
CSWAP[5, 5] = 0
CSWAP[5, 6] = 1
CSWAP[6, 5] = 1
CSWAP[6, 6] = 0
print(f"\nCSWAP|101⟩ = |110⟩: {np.allclose(CSWAP @ ket3(1, 0, 1), ket3(1, 1, 0))}") # True
print(f"CSWAP|011⟩ = |011⟩: {np.allclose(CSWAP @ ket3(0, 1, 1), ket3(0, 1, 1))}") # True
print(f"CSWAP is unitary: {is_unitary(CSWAP)}") # True
Operator Norm and Gate Distance
When you want to know how “close” two quantum gates are, you need a distance measure. The operator norm (also called the spectral norm) of the difference U - V gives a rigorous answer:
‖U - V‖ = largest singular value of (U - V)
Two gates are physically close if their operator norm distance is small. This metric matters for gate synthesis, where you approximate an ideal gate using a finite set of available gates, and for error analysis, where you quantify how much a noisy gate deviates from the ideal.
An important subtlety: two gates that differ only by a global phase are physically identical, even though their operator norm distance may be nonzero. The T gate and Rz(π/4) are a classic example.
# The T gate and Rz(pi/4) differ by a global phase
T_gate = np.array([[1, 0], [0, np.exp(1j * np.pi / 4)]], dtype=complex)
Rz_pi4 = Rz(np.pi / 4)
# Direct comparison: they are NOT equal
print(f"T == Rz(π/4): {np.allclose(T_gate, Rz_pi4)}") # False
# But they differ only by global phase e^{-iπ/8}
phase = np.exp(-1j * np.pi / 8)
print(f"T == e^{{-iπ/8}} Rz(π/4): {np.allclose(T_gate, phase * Rz_pi4)}") # True
# The operator norm of their difference after accounting for global phase is 0
diff = T_gate - phase * Rz_pi4
op_norm = np.linalg.norm(diff, ord=2) # largest singular value
print(f"‖T - phase*Rz(π/4)‖ = {op_norm:.2e}") # ~0
# Gate distance: Rx(π/2 + ε) vs Rx(π/2)
epsilon = 0.01
gate_ideal = Rx(np.pi / 2)
gate_noisy = Rx(np.pi / 2 + epsilon)
diff = gate_noisy - gate_ideal
op_norm = np.linalg.norm(diff, ord=2)
print(f"\n‖Rx(π/2 + {epsilon}) - Rx(π/2)‖ = {op_norm:.6f}")
print(f"ε/2 = {epsilon/2:.6f}")
# The operator norm is approximately ε/2 (to first order in ε),
# because the derivative of Rx(θ) with respect to θ has operator norm 1/2
For practical quantum computing, gate errors below about 10^{-3} to 10^{-4} are typically needed for error correction to work. The operator norm gives you a hardware-independent way to quantify whether you have met that threshold.
Building a Bell State with Matrix Math
A Bell state (maximally entangled two-qubit state) is created by applying H to qubit 0 then CNOT:
# Start with |00⟩
initial = np.kron(ket_0, ket_0) # [1, 0, 0, 0]
# Apply H to qubit 0: H ⊗ I
step1 = np.kron(H, I) @ initial
print(step1) # [0.70710678+0.j 0.+0.j 0.70710678+0.j 0.+0.j] = (|00⟩ + |10⟩) / √2
# Apply CNOT
bell_state = CNOT @ step1
print(bell_state) # [0.70710678+0.j 0.+0.j 0.+0.j 0.70710678+0.j] = (|00⟩ + |11⟩) / √2
The resulting state (|00⟩ + |11⟩) / √2 cannot be written as a tensor product of two individual qubit states; it is entangled. Both qubits are correlated: if you measure qubit 0 and get 0, qubit 1 is also 0; if you get 1, qubit 1 is also 1.
The Full Bell Basis Change Matrix
The Bell state circuit does something deeper than creating one entangled state. The full unitary matrix of the circuit (CNOT @ (H ⊗ I)) maps each computational basis state to a different Bell state:
|00⟩ → |Φ+⟩ = (|00⟩ + |11⟩) / √2
|01⟩ → |Ψ+⟩ = (|01⟩ + |10⟩) / √2
|10⟩ → |Φ-⟩ = (|00⟩ - |11⟩) / √2
|11⟩ → |Ψ-⟩ = (|01⟩ - |10⟩) / √2
In other words, this circuit is the Bell basis change matrix: it transforms the computational basis into the Bell basis.
# Build the full Bell circuit matrix
Bell_circuit = CNOT @ np.kron(H, I)
# Define the four computational basis inputs
state_00 = np.kron(ket_0, ket_0)
state_01 = np.kron(ket_0, ket_1)
state_10 = np.kron(ket_1, ket_0)
state_11 = np.kron(ket_1, ket_1)
# Apply the circuit to each
Phi_plus = Bell_circuit @ state_00
Psi_plus = Bell_circuit @ state_01
Phi_minus = Bell_circuit @ state_10
Psi_minus = Bell_circuit @ state_11
print("Bell states from computational basis inputs:")
print(f" |00⟩ → |Φ+⟩ = {np.round(Phi_plus, 3)}")
print(f" |01⟩ → |Ψ+⟩ = {np.round(Psi_plus, 3)}")
print(f" |10⟩ → |Φ-⟩ = {np.round(Phi_minus, 3)}")
print(f" |11⟩ → |Ψ-⟩ = {np.round(Psi_minus, 3)}")
# Verify the expected forms
inv_sqrt2 = 1 / np.sqrt(2)
assert np.allclose(Phi_plus, inv_sqrt2 * np.array([1, 0, 0, 1])) # (|00⟩ + |11⟩)/√2
assert np.allclose(Psi_plus, inv_sqrt2 * np.array([0, 1, 1, 0])) # (|01⟩ + |10⟩)/√2
assert np.allclose(Phi_minus, inv_sqrt2 * np.array([1, 0, 0, -1])) # (|00⟩ - |11⟩)/√2
assert np.allclose(Psi_minus, inv_sqrt2 * np.array([0, 1, -1, 0])) # (|01⟩ - |10⟩)/√2
print("\nAll four Bell states verified.")
# The inverse circuit (Bell measurement) is just the conjugate transpose
Bell_measure = Bell_circuit.conj().T
# It maps Bell states back to computational basis states
print("\nInverse (Bell measurement):")
print(f" |Φ+⟩ → {np.round(Bell_measure @ Phi_plus, 3)}") # [1, 0, 0, 0] = |00⟩
print(f" |Ψ+⟩ → {np.round(Bell_measure @ Psi_plus, 3)}") # [0, 1, 0, 0] = |01⟩
print(f" |Φ-⟩ → {np.round(Bell_measure @ Phi_minus, 3)}") # [0, 0, 1, 0] = |10⟩
print(f" |Ψ-⟩ → {np.round(Bell_measure @ Psi_minus, 3)}") # [0, 0, 0, 1] = |11⟩
This inverse operation is called a Bell measurement. It is the standard way to distinguish between the four Bell states: apply the inverse circuit, then measure in the computational basis. The outcome tells you which Bell state you had. Bell measurements are central to quantum teleportation and superdense coding.
Common Mistakes
When working with quantum gates as matrices, several pitfalls catch people repeatedly. Here are the most common ones and how to avoid them.
1. Wrong matrix multiplication order
For a circuit that applies gate A first and then gate B, the combined matrix is B @ A, not A @ B. The gate that acts first goes on the right, because matrix multiplication applies right-to-left to the state vector. This is counterintuitive since circuit diagrams read left-to-right.
# Circuit: H then X (read left to right)
# Matrix: X @ H (last gate on the left)
correct = X @ H @ ket_0
wrong = H @ X @ ket_0
print(f"Correct (X @ H)|0⟩: {np.round(correct, 3)}")
print(f"Wrong (H @ X)|0⟩: {np.round(wrong, 3)}")
print(f"Same? {np.allclose(correct, wrong)}") # False
2. Forgetting tensor products for multi-qubit systems
Applying H only to qubit 0 in a two-qubit system requires np.kron(H, I), not just H. A 2x2 matrix cannot multiply a 4-element state vector. This is a dimension error, and numpy will either raise an exception or silently produce a wrong answer via broadcasting.
two_qubit_state = np.kron(ket_0, ket_1) # |01⟩, a 4-element vector
# Correct: tensor product
H_on_q0 = np.kron(H, I) # 4x4 matrix
result = H_on_q0 @ two_qubit_state
print(f"Correct shape: {H_on_q0.shape} @ {two_qubit_state.shape} = {result.shape}")
# Wrong: using H directly would be a shape mismatch
# H @ two_qubit_state would raise ValueError: matmul: shapes (2,2) and (4,) not aligned
3. Confusing conjugate transpose with regular transpose
The adjoint (dagger) operation is the conjugate transpose: U† = U.conj().T. For real matrices, the transpose and conjugate transpose are the same. But quantum gates often have complex entries (Y, S, T all do), and using U.T instead of U.conj().T gives the wrong result.
# For the Y gate, transpose and conjugate transpose differ
print(f"Y.T:\n{Y.T}")
print(f"Y†:\n{Y.conj().T}")
print(f"Y.T == Y†: {np.allclose(Y.T, Y.conj().T)}") # False
# Only the conjugate transpose gives the correct inverse
print(f"Y† @ Y = I: {np.allclose(Y.conj().T @ Y, I)}") # True
print(f"Y.T @ Y = I: {np.allclose(Y.T @ Y, I)}") # False (gives -I)
4. Norm drift in long gate sequences
Each gate application introduces small floating-point errors. Over hundreds of gates, the state vector’s norm drifts away from 1.0. This is not a quantum effect; it is a numerical artifact. Periodically renormalizing the state keeps calculations accurate.
# Simulate norm drift over many random gates
state = ket_0.copy()
rng = np.random.default_rng(42)
for i in range(10000):
angle = rng.random() * 2 * np.pi
gate = Rx(angle)
state = gate @ state
print(f"Norm after 10000 gates: {np.linalg.norm(state):.15f}")
# May drift slightly from 1.0 due to floating point accumulation
# Fix: renormalize periodically
state = ket_0.copy()
for i in range(10000):
angle = rng.random() * 2 * np.pi
gate = Rx(angle)
state = gate @ state
if i % 1000 == 0:
state = state / np.linalg.norm(state)
print(f"Norm with periodic renormalization: {np.linalg.norm(state):.15f}")
5. Ignoring global phase in gate comparisons
Two gates that differ by a global phase are physically identical, but np.allclose will report them as different. For example, T² and S differ by a global phase:
# T² vs S
T_squared = T @ T
print(f"T² == S: {np.allclose(T_squared, S)}") # False!
# They differ by a global phase
# T = diag(1, e^{iπ/4}), so T² = diag(1, e^{iπ/2}) = diag(1, i) = S.
# Actually T² = S exactly. Let's check a trickier case.
# Rz(π/2) vs S: these differ by global phase
Rz_pi2 = Rz(np.pi / 2)
print(f"Rz(π/2):\n{np.round(Rz_pi2, 3)}")
print(f"S:\n{np.round(S, 3)}")
print(f"Rz(π/2) == S: {np.allclose(Rz_pi2, S)}") # False
# To compare up to global phase, check if U @ V† is proportional to I
def equal_up_to_phase(U, V, tol=1e-10):
"""Check if two unitaries are equal up to a global phase."""
product = U @ V.conj().T
# If U = e^{iφ} V, then U V† = e^{iφ} I
phase = product[0, 0]
return np.allclose(product, phase * np.eye(len(U)), atol=tol)
print(f"Rz(π/2) == S up to phase: {equal_up_to_phase(Rz_pi2, S)}") # True
print(f"T == Rz(π/4) up to phase: {equal_up_to_phase(T, Rz(np.pi/4))}") # True
Always use a phase-aware comparison when checking gate equivalences. Many textbooks and frameworks define gates with different global phase conventions, and naive equality checks will give false negatives.
What This Means for Quantum Computing
Every quantum circuit is, mathematically, a large unitary matrix acting on a quantum state vector. A circuit on n qubits operates on a 2^n-dimensional vector using a 2^n x 2^n unitary matrix. This exponential growth is both the promise and the challenge of quantum computing; representing the full matrix classically becomes infeasible around 50 qubits, but a quantum processor handles it physically at any scale.
Understanding gates as matrices makes it possible to reason about circuits analytically, compose gates on paper, and verify behavior before running on real hardware. The tools in this tutorial (eigenanalysis, gate decomposition, operator norms, projective measurement) form the mathematical foundation that all of quantum algorithm design builds on.
Was this tutorial helpful?