Complex Numbers in Quantum Computing
Why quantum computing uses complex numbers, how quantum amplitudes work, and the role of phase in quantum interference.
Quantum computing uses complex numbers everywhere. Quantum states, gate matrices, and measurement probabilities all involve complex arithmetic. This tutorial explains why complex numbers are necessary, how they work, and what role phase plays in quantum algorithms. Along the way you will see concrete Python code for every concept, so you can verify each claim yourself.
What Is a Complex Number?
A complex number has a real part and an imaginary part. The imaginary unit i is defined by the property that i squared equals negative one. A complex number z is written as:
z = a + bi
where a is the real part and b is the imaginary part.
import numpy as np
z1 = 2 + 3j # Python uses j for imaginary unit
z2 = 1 - 1j
print(z1 + z2) # (3+2j)
print(z1 * z2) # (5+1j)
print(z1.real, z1.imag) # 2.0, 3.0
The complex conjugate of z is written z* and flips the sign of the imaginary part:
print(z1.conjugate()) # (2-3j)
The magnitude (or absolute value) of a complex number measures its distance from zero:
print(abs(z1)) # 3.606 = sqrt(2^2 + 3^2)
The Polar Form and Euler’s Formula
Any complex number can be written in polar form using a magnitude r and an angle theta:
z = r * e^(i * theta) = r * (cos(theta) + i * sin(theta))
This is Euler’s formula, and it is one of the most useful identities in quantum computing. It connects complex exponentials to rotations.
r = 1.0
theta = np.pi / 4 # 45 degrees
z = r * np.exp(1j * theta)
print(z) # (0.707+0.707j)
print(abs(z)) # 1.0 -- magnitude unchanged
# The angle (phase) of z
print(np.angle(z)) # 0.785 rad = pi/4
Complex Multiplication Is Rotation and Scaling
In polar form, multiplication becomes transparent. Given z1 = r1 * e^(i * theta1) and z2 = r2 * e^(i * theta2):
z1 * z2 = r1 * r2 * e^(i * (theta1 + theta2))
The magnitudes multiply. The angles add. Multiplying by a unit-magnitude complex number e^(i * phi) is a pure rotation by angle phi, with no change in magnitude.
# Demonstrate: multiplication = rotation + scaling
z1 = 2 * np.exp(1j * np.pi / 6) # magnitude 2, angle 30 degrees
z2 = 3 * np.exp(1j * np.pi / 3) # magnitude 3, angle 60 degrees
product = z1 * z2
print(f"Magnitudes: {abs(z1):.1f} * {abs(z2):.1f} = {abs(product):.1f}")
# Magnitudes: 2.0 * 3.0 = 6.0
print(f"Angles: {np.degrees(np.angle(z1)):.0f} + {np.degrees(np.angle(z2)):.0f} "
f"= {np.degrees(np.angle(product)):.0f} degrees")
# Angles: 30 + 60 = 90 degrees
This is exactly what phase gates do in quantum computing. The S gate multiplies the |1> amplitude by e^(i * pi/2) = i, which rotates that amplitude by 90 degrees in the complex plane without changing its magnitude:
S = np.array([[1, 0],
[0, 1j]], dtype=complex)
# Start with |+> = (|0> + |1>) / sqrt(2)
psi = np.array([1/np.sqrt(2), 1/np.sqrt(2)], dtype=complex)
after_S = S @ psi
print(f"|0> amplitude: {after_S[0]:.4f}") # 0.7071+0.0000j (unchanged)
print(f"|1> amplitude: {after_S[1]:.4f}") # 0.0000+0.7071j (rotated 90 degrees)
print(f"|1> phase: {np.degrees(np.angle(after_S[1])):.0f} degrees") # 90 degrees
print(f"|1> magnitude: {abs(after_S[1]):.4f}") # 0.7071 (unchanged)
The probabilities remain 50/50, but the phase of the |1> component has shifted. That phase shift will become visible when a subsequent gate converts it into a probability difference.
Why Real Numbers Are Not Enough
You might wonder whether complex numbers are truly necessary. Some quantum gates are purely real. The Ry rotation gate, for example, has only real entries:
Ry(theta) = [[cos(theta/2), -sin(theta/2)],
[sin(theta/2), cos(theta/2)]]
def Ry(theta):
c = np.cos(theta / 2)
s = np.sin(theta / 2)
return np.array([[c, -s],
[s, c]], dtype=complex)
print(Ry(np.pi / 2))
# [[ 0.707 -0.707]
# [ 0.707 0.707]]
This gate is a real orthogonal matrix (an element of SO(2)). It can rotate a qubit’s amplitudes, and the entries are all real. So why do we need complex numbers at all?
The Y gate requires complex entries
The Pauli Y gate performs a specific rotation that cannot be expressed with real numbers alone:
Y = np.array([[0, -1j],
[1j, 0]], dtype=complex)
# Verify it is unitary: Y† Y = I
print(Y.conj().T @ Y)
# [[1.+0.j 0.+0.j]
# [0.+0.j 1.+0.j]]
If you try to replace the i entries with real numbers while keeping the gate unitary, you get a different gate with different physics.
The full group of single-qubit unitaries requires three parameters
The deeper reason is mathematical. The most general single-qubit unitary gate belongs to the group SU(2) and requires three real parameters (theta, phi, lambda):
U(theta, phi, lambda) = [[cos(theta/2), -e^(i*lambda) * sin(theta/2)],
[e^(i*phi) * sin(theta/2), e^(i*(phi+lambda)) * cos(theta/2)]]
def U_gate(theta, phi, lam):
"""General single-qubit unitary in SU(2)."""
return np.array([
[np.cos(theta/2), -np.exp(1j * lam) * np.sin(theta/2)],
[np.exp(1j * phi) * np.sin(theta/2), np.exp(1j * (phi + lam)) * np.cos(theta/2)]
], dtype=complex)
# Example: theta=pi/2, phi=pi/4, lambda=pi/3
U = U_gate(np.pi/2, np.pi/4, np.pi/3)
# Verify unitarity
identity_check = U.conj().T @ U
print(np.allclose(identity_check, np.eye(2))) # True
Real orthogonal 2x2 matrices (SO(2)) have only one parameter: the rotation angle. That gives you a one-dimensional family of gates. SU(2) is three-dimensional. Complex numbers give you access to the full space of single-qubit operations that quantum mechanics allows.
Without complex numbers, you cannot reach arbitrary points on the Bloch sphere from an arbitrary starting state. You would be restricted to rotations in a single plane.
Why Quantum Computing Needs Complex Numbers
A quantum state is a vector of complex amplitudes. For a qubit:
|psi> = alpha |0> + beta |1>
where alpha and beta are complex numbers satisfying |alpha|^2 + |beta|^2 = 1.
The probabilities of measuring each outcome are the squared magnitudes of the amplitudes. The imaginary parts do not appear in the measurement probabilities directly, but they are essential for how quantum states evolve and interfere.
# A state with complex amplitudes
alpha = 1 / np.sqrt(2)
beta = 1j / np.sqrt(2) # imaginary amplitude
prob_0 = abs(alpha)**2
prob_1 = abs(beta)**2
print(f"P(0) = {prob_0:.3f}") # 0.500
print(f"P(1) = {prob_1:.3f}") # 0.500
print(f"Sum = {prob_0 + prob_1:.3f}") # 1.000
Even though the probabilities are 50/50 just like the |+> state, this state is different. Its phase structure will produce different results when gates are applied.
The Bloch Sphere and Complex Amplitudes
Every single-qubit pure state can be written as:
|psi> = cos(theta/2) |0> + e^(i*phi) * sin(theta/2) |1>
where theta and phi are the polar and azimuthal angles on the Bloch sphere. The global phase has been factored out, leaving two real parameters from the complex amplitude ratio.
The key mapping is:
- theta = 0 (north pole): the state is |0>. The azimuthal angle phi is undefined because sin(0) = 0, so the e^(i*phi) factor vanishes entirely.
- theta = pi (south pole): the state is |1>.
- theta = pi/2, phi = 0: the state is |+> = (|0> + |1>) / sqrt(2).
- theta = pi/2, phi = pi: the state is |-> = (|0> - |1>) / sqrt(2).
- theta = pi/2, phi = pi/2: the state is |+i> = (|0> + i|1>) / sqrt(2).
Every point on the Bloch sphere corresponds to a unique complex amplitude ratio beta/alpha = e^(i*phi) * tan(theta/2).
def bloch_to_statevector(theta, phi):
"""Convert Bloch sphere coordinates (theta, phi) to a state vector."""
alpha = np.cos(theta / 2)
beta = np.exp(1j * phi) * np.sin(theta / 2)
return np.array([alpha, beta], dtype=complex)
# North pole: |0>
print("North pole:", bloch_to_statevector(0, 0))
# [1.+0.j 0.+0.j]
# South pole: |1>
print("South pole:", bloch_to_statevector(np.pi, 0))
# [0. 1.+0.j] (within numerical precision)
# Equator, phi=0: |+>
psi_plus = bloch_to_statevector(np.pi/2, 0)
print("|+>:", psi_plus)
# [0.707+0.j 0.707+0.j]
# Equator, phi=pi/2: |+i>
psi_plus_i = bloch_to_statevector(np.pi/2, np.pi/2)
print("|+i>:", psi_plus_i)
# [0.707+0.j 0.+0.707j]
# Verify normalization for all states
for name, (th, ph) in [("pole", (0.3, 1.7)), ("equator", (np.pi/2, 2.1))]:
psi = bloch_to_statevector(th, ph)
print(f"{name}: |alpha|^2 + |beta|^2 = {np.abs(psi[0])**2 + np.abs(psi[1])**2:.6f}")
# Always 1.000000
Inner Products and Probability from Complex Amplitudes
The inner product of two quantum states |psi> and |phi> is:
<psi|phi> = sum_i psi_i* phi_i
where psi_i* denotes the complex conjugate of psi_i. The conjugate on the left vector (the “bra”) is not optional. It is what makes the inner product positive-definite, meaning <psi|psi> is always a non-negative real number.
# Define two states
psi = np.array([1/np.sqrt(2), 1j/np.sqrt(2)], dtype=complex)
phi = np.array([1/np.sqrt(2), 1/np.sqrt(2)], dtype=complex)
# Inner product: <psi|phi> = psi† . phi
inner = np.vdot(psi, phi) # np.vdot conjugates the first argument
print(f"<psi|phi> = {inner:.4f}") # (0.5000-0.5000j)
# Normalization check: <psi|psi> = 1
norm = np.vdot(psi, psi)
print(f"<psi|psi> = {norm:.4f}") # (1.0000+0.0000j)
Why the conjugate matters
Without the conjugate, the “inner product” of a state with itself could be zero or even negative for complex amplitudes. Consider psi = [1, i] / sqrt(2):
psi = np.array([1, 1j], dtype=complex) / np.sqrt(2)
# Wrong: no conjugate (just a dot product)
wrong_norm = psi @ psi
print(f"Without conjugate: {wrong_norm:.4f}") # (0.0000+0.0000j) -- zero!
# Correct: with conjugate
correct_norm = psi.conj() @ psi
print(f"With conjugate: {correct_norm:.4f}") # (1.0000+0.0000j) -- one, as expected
The “wrong” version gives zero because 11 + ii = 1 + (-1) = 0. The conjugate fixes this: 11 + (-i)(i) = 1 + 1 = 2, and after normalization, we get 1.
The Born rule from inner products
The probability of measuring a state |psi> in the computational basis state |0> is:
P(0) = |<0|psi>|^2 = |alpha*|^2 = |alpha|^2
where alpha is the |0> component of |psi>. The conjugate on alpha and the subsequent squaring of the magnitude produce a real, non-negative probability.
psi = np.array([0.6 + 0.2j, 0.3 + 0.7j], dtype=complex)
# Normalize
psi = psi / np.linalg.norm(psi)
# Basis states
ket_0 = np.array([1, 0], dtype=complex)
ket_1 = np.array([0, 1], dtype=complex)
# Probabilities via inner products
p0 = np.abs(np.vdot(ket_0, psi))**2
p1 = np.abs(np.vdot(ket_1, psi))**2
print(f"P(0) = {p0:.4f}")
print(f"P(1) = {p1:.4f}")
print(f"Sum = {p0 + p1:.4f}") # 1.0000
Phase: What It Is and Why It Matters
The phase of a complex amplitude is its angle in the complex plane. Two states can have the same measurement probabilities but different phases, and those phases determine how the states respond to gates.
There are two kinds of phase to distinguish:
Global phase is an overall factor e^(i*theta) applied to the entire state. It has no observable effect and can always be ignored:
psi1 = np.array([1/np.sqrt(2), 1/np.sqrt(2)], dtype=complex)
psi2 = np.exp(1j * np.pi/4) * psi1 # global phase applied
# Probabilities are identical
print(abs(psi1)**2) # [0.5, 0.5]
print(abs(psi2)**2) # [0.5, 0.5]
Relative phase is a phase difference between components of a superposition. Relative phase is physically meaningful because it affects interference.
# |+> state: both amplitudes in phase
psi_plus = np.array([1/np.sqrt(2), 1/np.sqrt(2)], dtype=complex)
# |-> state: amplitudes out of phase by pi
psi_minus = np.array([1/np.sqrt(2), -1/np.sqrt(2)], dtype=complex)
# Apply Hadamard to both
H = (1/np.sqrt(2)) * np.array([[1, 1], [1, -1]], dtype=complex)
print(H @ psi_plus) # [1.+0.j, 0.+0.j] -- collapses to |0>
print(H @ psi_minus) # [0.+0.j, 1.+0.j] -- collapses to |1>
The Hadamard gate converts relative phase into a difference in measurement probabilities. This is the essence of quantum interference.
Important subtlety: Global phase is genuinely unobservable for isolated qubits. However, relative phase between qubits in an entangled state IS observable. When qubit A is entangled with qubit B, a “global” phase on qubit A alone becomes a relative phase of the joint system and can produce measurable effects.
Interference: How Phase Affects Computation
Quantum interference happens when probability amplitudes add together. Amplitudes with the same phase add constructively (making the outcome more likely). Amplitudes with opposite phases cancel destructively (making the outcome less likely).
# Constructive interference
a = 1/np.sqrt(2)
b = 1/np.sqrt(2)
combined = a + b # amplitudes add up
print(abs(combined)**2) # 2.0 -- but this would need renormalization
# Destructive interference
c = 1/np.sqrt(2)
d = -1/np.sqrt(2)
cancelled = c + d
print(abs(cancelled)**2) # 0.0 -- complete cancellation
In a real quantum algorithm, the circuit is designed so that amplitudes leading to wrong answers cancel and amplitudes leading to correct answers reinforce. This is how Grover’s algorithm and quantum Fourier transform-based algorithms achieve their speedups.
Interference in the Double-Slit Analogy
In a double-slit experiment, a particle can take two paths to reach a detector. Each path contributes a complex amplitude. The total amplitude at the detector is the sum of both path amplitudes, and the probability is the squared magnitude of that sum.
This is the same principle at work in quantum circuits. Each computational path through a circuit contributes a complex amplitude to a given output. The amplitudes from all paths sum before squaring to get a probability.
The Deutsch algorithm provides a clean example. Consider a function f that maps one bit to one bit. The algorithm determines whether f is constant (f(0) = f(1)) or balanced (f(0) != f(1)) with a single query:
H = (1/np.sqrt(2)) * np.array([[1, 1], [1, -1]], dtype=complex)
# After the Deutsch circuit, the amplitude of measuring |0> on the first qubit is:
# ((-1)^f(0) + (-1)^f(1)) / 2
# Case 1: constant function, f(0)=0, f(1)=0
f0, f1 = 0, 0
amp_0 = ((-1)**f0 + (-1)**f1) / 2
amp_1 = ((-1)**f0 - (-1)**f1) / 2
print(f"Constant: amp(|0>)={amp_0:.2f}, amp(|1>)={amp_1:.2f}")
print(f" P(|0>)={abs(amp_0)**2:.2f}, P(|1>)={abs(amp_1)**2:.2f}")
# Constant: amp(|0>)=1.00, amp(|1>)=0.00
# Constructive interference on |0>, destructive on |1>
# Case 2: balanced function, f(0)=0, f(1)=1
f0, f1 = 0, 1
amp_0 = ((-1)**f0 + (-1)**f1) / 2
amp_1 = ((-1)**f0 - (-1)**f1) / 2
print(f"Balanced: amp(|0>)={amp_0:.2f}, amp(|1>)={amp_1:.2f}")
print(f" P(|0>)={abs(amp_0)**2:.2f}, P(|1>)={abs(amp_1)**2:.2f}")
# Balanced: amp(|0>)=0.00, amp(|1>)=1.00
# Destructive interference on |0>, constructive on |1>
For the constant function, both paths contribute the same sign, so they add constructively. For the balanced function, the paths contribute opposite signs and cancel. Interference converts a phase difference into a deterministic measurement outcome.
Complex Numbers in Quantum Gates
All standard quantum gates have complex entries. The Y gate and the S gate illustrate this:
Y = np.array([[0, -1j],
[1j, 0]], dtype=complex)
S = np.array([[1, 0],
[0, 1j]], dtype=complex)
# S gate adds a phase of i to the |1> component
psi = np.array([1/np.sqrt(2), 1/np.sqrt(2)], dtype=complex)
after_S = S @ psi
print(after_S) # [0.707+0.j, 0.+0.707j]
# The |1> component now has a phase of i = e^(i*pi/2)
The T gate (pi/8 gate) adds a phase of e^(i*pi/4) to the |1> component:
T = np.array([[1, 0],
[0, np.exp(1j * np.pi / 4)]], dtype=complex)
after_T = T @ psi
print(after_T) # [0.707+0.j, 0.5+0.5j]
print(np.angle(after_T[1])) # 0.785 rad = pi/4
Phase Kickback
Phase kickback is one of the most important phenomena in quantum algorithms. It occurs when applying a controlled-U gate and the target qubit is an eigenstate of U. Instead of the target qubit changing, the control qubit acquires a phase.
Here is the setup: U is a unitary with eigenstate |u> and eigenvalue e^(iphi). That is, U|u> = e^(iphi)|u>. When we apply controlled-U with the control qubit in superposition and the target in |u>:
CU (|0> + |1>)/sqrt(2) |u> = (|0>|u> + e^(i*phi)|1>|u>) / sqrt(2)
= ((|0> + e^(i*phi)|1>) / sqrt(2)) |u>
The phase “kicks back” onto the control qubit. The target qubit is unchanged.
# Concrete example: CZ gate, target in |1> (eigenstate of Z with eigenvalue -1)
# Z gate
Z = np.array([[1, 0], [0, -1]], dtype=complex)
# |1> is an eigenstate of Z with eigenvalue -1 = e^(i*pi)
eigenstate = np.array([0, 1], dtype=complex)
print(f"Z|1> = {Z @ eigenstate}") # [0, -1] = -1 * |1> = e^(i*pi) * |1>
# Control qubit starts in |+>
control = np.array([1/np.sqrt(2), 1/np.sqrt(2)], dtype=complex)
# CZ gate as a 4x4 matrix
# CZ = |0><0| x I + |1><1| x Z
I = np.eye(2, dtype=complex)
proj_0 = np.array([[1, 0], [0, 0]], dtype=complex) # |0><0|
proj_1 = np.array([[0, 0], [0, 1]], dtype=complex) # |1><1|
CZ = np.kron(proj_0, I) + np.kron(proj_1, Z)
# Full initial state: control x target
initial = np.kron(control, eigenstate)
print(f"Initial state: {initial}")
# [0. 0.707 0. 0.707]
# Apply CZ
final = CZ @ initial
print(f"Final state: {final}")
# [0. 0.707 0. -0.707]
# Factor out the target (it's still |1>):
# The control qubit is now (|0> - |1>)/sqrt(2) = |->
# Phase e^(i*pi) = -1 has kicked back to the control qubit
control_after = np.array([final[1], final[3]]) # extract control amplitudes
print(f"Control after: {control_after}")
# [0.707, -0.707] = |->
# Verify: this is |-> up to normalization
minus_state = np.array([1/np.sqrt(2), -1/np.sqrt(2)], dtype=complex)
print(f"Is |->? {np.allclose(control_after, minus_state)}") # True
The control qubit started in |+> and ended in |->. The Z eigenvalue of -1 = e^(i*pi) shifted the relative phase of the control qubit by pi. This is the mechanism behind quantum phase estimation and many other quantum algorithms.
Complex Exponentials in Quantum Algorithms
The quantum Fourier transform (QFT) is built from roots of unity. The Nth roots of unity are the N complex numbers:
omega^k = e^(2*pi*i*k/N) for k = 0, 1, ..., N-1
These N points are equally spaced on the unit circle. Each has magnitude 1 (pure phase, no amplitude change) and a distinct angle.
# Roots of unity for N=8
N = 8
roots = [np.exp(2j * np.pi * k / N) for k in range(N)]
print("N=8 roots of unity:")
for k, w in enumerate(roots):
print(f" omega^{k} = {w:.4f} "
f"(magnitude={abs(w):.4f}, angle={np.degrees(np.angle(w)):7.1f} deg)")
The DFT matrix has entries F[j,k] = omega^(jk) / sqrt(N). For N=4:
N = 4
omega = np.exp(2j * np.pi / N) # = i
# Build the 4x4 DFT matrix
F = np.zeros((N, N), dtype=complex)
for j in range(N):
for k in range(N):
F[j, k] = omega**(j * k) / np.sqrt(N)
print("DFT matrix (N=4):")
print(np.round(F, 3))
# [[ 0.5+0.j 0.5+0.j 0.5+0.j 0.5+0.j ]
# [ 0.5+0.j 0. +0.5j -0.5+0.j 0. -0.5j]
# [ 0.5+0.j -0.5+0.j 0.5+0.j -0.5+0.j ]
# [ 0.5+0.j 0. -0.5j -0.5+0.j 0. +0.5j]]
# Verify unitarity: F† F = I
print("\nF† F = I?", np.allclose(F.conj().T @ F, np.eye(N))) # True
# Verify all entries have magnitude 1/sqrt(N)
print("All entries magnitude 1/sqrt(N)?",
np.allclose(np.abs(F), 1/np.sqrt(N))) # True
Because every entry of the DFT matrix has the same magnitude (1/sqrt(N)), the QFT preserves normalization. The complex phases are what encode the frequency information. This is the mathematical foundation of Shor’s factoring algorithm and quantum phase estimation.
Measuring Complex Amplitudes: Quantum State Tomography
You cannot directly measure the complex amplitudes of a quantum state. A measurement in the computational basis yields only probabilities, which are the squared magnitudes |alpha|^2 and |beta|^2. All phase information is lost.
To reconstruct the full complex state, you need to measure in multiple bases. This procedure is called quantum state tomography.
For a single qubit |psi> = alpha|0> + beta|1>, three sets of measurements suffice:
Step 1: Measure in the Z basis to get |alpha|^2 and |beta|^2.
Step 2: Apply H, then measure in Z to get information about Re(alpha* beta). After H, the state becomes [(alpha + beta), (alpha - beta)] / sqrt(2). The probability of outcome 0 is:
P_X(0) = |alpha + beta|^2 / 2 = (|alpha|^2 + |beta|^2 + 2 Re(alpha* beta)) / 2
So Re(alpha* beta) = P_X(0) - 1/2.
Step 3: Apply S-dagger then H, then measure in Z to get Im(alpha* beta). S-dagger maps |1> to -i|1>, so after S-dagger, beta becomes -i*beta. Then H gives the probability:
P_Y(0) = |alpha + (-i*beta)|^2 / 2 = (|alpha|^2 + |beta|^2 + 2 Re(alpha* (-i*beta))) / 2
= (1 + 2 Im(alpha* beta)) / 2
So Im(alpha* beta) = P_Y(0) - 1/2.
def tomography_reconstruct(pz0, px0, py0):
"""
Reconstruct a single-qubit state from tomography probabilities.
pz0: probability of |0> in Z measurement
px0: probability of |0> after H (X measurement)
py0: probability of |0> after S†H (Y measurement)
Returns the reconstructed state vector (up to global phase).
"""
# From Z measurement
abs_alpha = np.sqrt(pz0)
abs_beta = np.sqrt(1 - pz0)
# Choose alpha real and positive (fixing global phase)
alpha = abs_alpha
# From X and Y measurements
re_part = px0 - 0.5 # Re(alpha* beta)
im_part = py0 - 0.5 # Im(alpha* beta)
if abs_alpha > 1e-10:
# alpha* beta = alpha * beta (since alpha is real)
# So beta = (re_part + i*im_part) / alpha
beta = (re_part + 1j * im_part) / alpha
else:
# alpha is ~0, state is ~|1>, phase of beta is global phase
beta = abs_beta
state = np.array([alpha, beta], dtype=complex)
state = state / np.linalg.norm(state) # renormalize
return state
# Test: reconstruct the state |psi> = cos(pi/8)|0> + e^(i*pi/3)*sin(pi/8)|1>
theta_true = np.pi / 4
phi_true = np.pi / 3
alpha_true = np.cos(theta_true / 2)
beta_true = np.exp(1j * phi_true) * np.sin(theta_true / 2)
psi_true = np.array([alpha_true, beta_true], dtype=complex)
print(f"True state: alpha={alpha_true:.4f}, beta={beta_true:.4f}")
# Simulate ideal tomography measurements
H = (1/np.sqrt(2)) * np.array([[1, 1], [1, -1]], dtype=complex)
Sdg = np.array([[1, 0], [0, -1j]], dtype=complex)
# Z measurement
pz0 = np.abs(psi_true[0])**2
# X measurement (apply H first)
psi_x = H @ psi_true
px0 = np.abs(psi_x[0])**2
# Y measurement (apply S†H first)
psi_y = H @ Sdg @ psi_true
py0 = np.abs(psi_y[0])**2
print(f"Tomography probabilities: pz0={pz0:.4f}, px0={px0:.4f}, py0={py0:.4f}")
# Reconstruct
psi_recon = tomography_reconstruct(pz0, px0, py0)
print(f"Reconstructed: alpha={psi_recon[0]:.4f}, beta={psi_recon[1]:.4f}")
# Check: states should match up to global phase
overlap = np.abs(np.vdot(psi_true, psi_recon))**2
print(f"Fidelity: {overlap:.6f}") # 1.000000
This demonstrates that complex amplitudes, while not directly observable, leave measurable signatures when you probe the state from different angles.
Complex Numbers in Density Matrices
For mixed states (noisy qubits or subsystems of entangled states), the density matrix provides a complete description:
rho = sum_i p_i |psi_i><psi_i|
A density matrix is Hermitian (rho = rho†), meaning it equals its own conjugate transpose. This forces the diagonal entries to be real and the off-diagonal entries to be complex conjugates of each other.
# Pure state |+> = (|0> + |1>) / sqrt(2)
psi_plus = np.array([[1/np.sqrt(2)], [1/np.sqrt(2)]], dtype=complex)
rho_plus = psi_plus @ psi_plus.conj().T
print("Density matrix for |+>:")
print(rho_plus)
# [[0.5 0.5]
# [0.5 0.5]]
# Off-diagonal entries are 0.5: full coherence
The off-diagonal entry rho_01 = <0|rho|1> encodes the coherence between |0> and |1>. Its magnitude tells you how much quantum superposition remains, and its phase tells you which superposition.
# Pure state |+i> = (|0> + i|1>) / sqrt(2)
psi_plus_i = np.array([[1/np.sqrt(2)], [1j/np.sqrt(2)]], dtype=complex)
rho_plus_i = psi_plus_i @ psi_plus_i.conj().T
print("Density matrix for |+i>:")
print(rho_plus_i)
# [[0.5 -0.5j]
# [0.5j 0.5 ]]
# Off-diagonal entries are imaginary: the phase is pi/2
# Completely mixed state (maximally noisy)
rho_mixed = np.array([[0.5, 0], [0, 0.5]], dtype=complex)
print("\nCompletely mixed state:")
print(rho_mixed)
# [[0.5 0. ]
# [0. 0.5]]
# Off-diagonal entries are zero: no coherence
Decoherence is precisely the process of the off-diagonal entries decaying toward zero. A qubit that starts in a coherent superposition gradually loses its complex off-diagonal entries, and with them, its ability to interfere. This is why quantum computers need error correction: without it, the complex structure that makes quantum computation powerful degrades over time.
# Simulating dephasing: off-diagonal entries decay exponentially
def dephased_state(rho, gamma):
"""Apply dephasing channel with parameter gamma in [0, 1].
gamma=0: no dephasing, gamma=1: complete dephasing."""
result = rho.copy()
result[0, 1] *= (1 - gamma)
result[1, 0] *= (1 - gamma)
return result
rho_initial = rho_plus.copy()
for gamma in [0, 0.25, 0.5, 0.75, 1.0]:
rho_t = dephased_state(rho_initial, gamma)
print(f"gamma={gamma:.2f}: rho_01={rho_t[0,1]:.3f}, "
f"coherence={abs(rho_t[0,1]):.3f}")
# gamma=0.00: rho_01=0.500, coherence=0.500
# gamma=0.25: rho_01=0.375, coherence=0.375
# gamma=0.50: rho_01=0.250, coherence=0.250
# gamma=0.75: rho_01=0.125, coherence=0.125
# gamma=1.00: rho_01=0.000, coherence=0.000
Worked Example: Phase Estimation Intuition
Quantum phase estimation (QPE) extracts the eigenvalue phase of a unitary operator. Let us walk through a minimal example with a single ancilla qubit.
Suppose U has eigenvalue e^(ipi/2) = i (eigenphase 1/4, since pi/2 = 2pi * 1/4). The ancilla starts in |+>, and the target is in the eigenstate of U.
After controlled-U, the ancilla state is:
(|0> + e^(i*pi/2)|1>) / sqrt(2) = (|0> + i|1>) / sqrt(2)
Now we apply H (which acts as a 1-qubit inverse QFT) to the ancilla:
# Ancilla state after controlled-U
ancilla = np.array([1/np.sqrt(2), 1j/np.sqrt(2)], dtype=complex)
# Apply H (inverse QFT for 1 qubit)
H = (1/np.sqrt(2)) * np.array([[1, 1], [1, -1]], dtype=complex)
measured = H @ ancilla
print(f"After H: {measured}")
# After H: [(0.5+0.5j), (0.5-0.5j)]
# Measurement probabilities
p0 = np.abs(measured[0])**2
p1 = np.abs(measured[1])**2
print(f"P(0) = {p0:.4f}") # 0.5000
print(f"P(1) = {p1:.4f}") # 0.5000
With a single ancilla qubit, both outcomes have equal probability. The single ancilla can only distinguish eigenphases of 0 and 1/2 perfectly. The eigenphase 1/4 falls between these, so neither outcome is certain.
With more ancilla qubits, the resolution improves. Using n ancilla qubits gives 2^n bins for the phase, allowing you to resolve eigenphases to n bits of precision:
def simulate_qpe(eigenphase, n_ancilla):
"""
Simulate ideal QPE outcome probabilities for a given eigenphase.
eigenphase: the phase phi where eigenvalue = e^(2*pi*i*phi), in [0, 1)
n_ancilla: number of ancilla qubits
"""
N = 2**n_ancilla
# After controlled-U operations and inverse QFT,
# the probability of measuring outcome m is:
# P(m) = |sum_{k=0}^{N-1} e^(2*pi*i*k*(phi - m/N)) / N|^2
probs = np.zeros(N)
for m in range(N):
amp = sum(np.exp(2j * np.pi * k * (eigenphase - m/N))
for k in range(N)) / N
probs[m] = np.abs(amp)**2
return probs
eigenphase = 0.25 # eigenvalue = e^(i*pi/2) = i
for n in [1, 2, 3, 4]:
probs = simulate_qpe(eigenphase, n)
N = 2**n
best = np.argmax(probs)
print(f"n={n} ancilla ({N:2d} bins): "
f"best outcome = {best}/{N} = {best/N:.4f}, "
f"P = {probs[best]:.4f}")
# n=1 ancilla ( 2 bins): best outcome = 0/2 = 0.0000, P = 0.5000
# n=2 ancilla ( 4 bins): best outcome = 1/4 = 0.2500, P = 1.0000
# n=3 ancilla ( 8 bins): best outcome = 2/8 = 0.2500, P = 1.0000
# n=4 ancilla (16 bins): best outcome = 4/16 = 0.2500, P = 1.0000
With 2 or more ancilla qubits, the eigenphase 1/4 is resolved exactly (because 1/4 is a multiple of 1/N for N >= 4). For eigenphases that do not align with a bin, more ancilla qubits provide better approximations, with the probability concentrating around the nearest bin.
NumPy Reference for Complex Arithmetic in Quantum Computing
Here is a compact reference of the complex number operations you will use repeatedly in quantum computing code:
import numpy as np
z = 0.5 + 0.5j
# --- Scalar operations ---
np.conj(z) # Complex conjugate: (0.5-0.5j)
np.abs(z) # Magnitude: 0.7071
np.angle(z) # Phase angle in radians: 0.7854 (pi/4)
np.exp(1j * np.pi) # Euler's formula: e^(i*pi) = -1
# --- Array operations ---
psi = np.array([1/np.sqrt(2), 1j/np.sqrt(2)], dtype=complex)
phi = np.array([1/np.sqrt(2), 1/np.sqrt(2)], dtype=complex)
np.conj(psi) # Element-wise conjugate
np.abs(psi)**2 # Probabilities: [0.5, 0.5]
np.angle(psi) # Phase of each amplitude
# --- Inner product (conjugates the first argument) ---
np.vdot(psi, phi) # <psi|phi> = psi† . phi
# --- Matrix-vector product ---
H = (1/np.sqrt(2)) * np.array([[1, 1], [1, -1]], dtype=complex)
H @ psi # Apply gate H to state psi
# --- Hermitian conjugate (dagger) ---
H.conj().T # H† = conjugate transpose
# --- Outer product (for density matrices) ---
rho = np.outer(psi, np.conj(psi)) # |psi><psi|
# --- Expectation value <psi|A|psi> ---
A = np.array([[1, 0], [0, -1]], dtype=complex) # Pauli Z
expect = psi.conj() @ A @ psi # Note: conj() on the left vector
print(f"<Z> = {expect.real:.4f}")
# --- Verify unitarity ---
print(np.allclose(H.conj().T @ H, np.eye(2))) # True
Common Mistakes with Complex Numbers in Quantum Code
1. Forgetting np.abs() for arrays
Python’s built-in abs() works on individual complex numbers but np.abs() is preferred for arrays because it handles element-wise operations cleanly:
psi = np.array([0.5 + 0.5j, 0.5 - 0.5j])
# Both work for arrays, but numpy is idiomatic
print(abs(psi)) # [0.7071, 0.7071] -- works, calls np.abs internally
print(np.abs(psi)) # [0.7071, 0.7071] -- explicit and preferred
2. Confusing conjugate methods
Both .conjugate() (method) and np.conj() (function) compute the complex conjugate. For arrays, the numpy function reads more clearly:
z = 3 + 4j
print(z.conjugate()) # (3-4j)
print(np.conj(z)) # (3-4j)
arr = np.array([1+2j, 3-4j])
print(np.conj(arr)) # [1-2j, 3+4j] -- cleaner for arrays
3. Misusing Python exponentiation for Euler’s formula
np.exp(1j * theta) correctly computes e^(itheta) using Euler’s formula. A common mistake is writing np.e ** (1j * theta), which computes (2.71828…)^(itheta) using Python’s power operator. While mathematically equivalent, the np.exp form is numerically safer and universally understood:
theta = np.pi / 4
# Correct and preferred
result1 = np.exp(1j * theta)
# Also correct but avoid for clarity
result2 = np.e ** (1j * theta)
print(np.allclose(result1, result2)) # True, but use np.exp
4. Forgetting the conjugate in expectation values
The expectation value <psi|A|psi> requires conjugating psi on the left side. Omitting the conjugate produces wrong results for states with complex amplitudes:
A = np.array([[1, 0], [0, -1]], dtype=complex) # Pauli Z
psi = np.array([1/np.sqrt(2), 1j/np.sqrt(2)], dtype=complex)
# Wrong: no conjugate
wrong = psi @ A @ psi
print(f"Wrong: {wrong:.4f}") # (0.0000+0.0000j) -- not even real!
# Correct: conjugate on the left
correct = psi.conj() @ A @ psi
print(f"Correct: {correct:.4f}") # (0.0000+0.0000j) -- happens to be 0 here too
# More dramatic example where it matters
psi2 = np.array([np.sqrt(0.8), 1j * np.sqrt(0.2)], dtype=complex)
wrong2 = psi2 @ A @ psi2
correct2 = psi2.conj() @ A @ psi2
print(f"Wrong: {wrong2:.4f}") # (0.6000+0.0000j) -- wrong
print(f"Correct: {correct2:.4f}") # (0.6000+0.0000j) -- same here because Z is diagonal
# But for non-diagonal operators the difference is stark:
X = np.array([[0, 1], [1, 0]], dtype=complex)
wrong3 = psi2 @ X @ psi2
correct3 = psi2.conj() @ X @ psi2
print(f"<X> wrong: {wrong3:.4f}") # May have imaginary part
print(f"<X> correct: {correct3:.4f}") # Always real for Hermitian operators
The rule is simple: expectation values of Hermitian operators are always real. If your result has an imaginary part, check that you conjugated the left vector.
5. Assuming global phase is always irrelevant
Global phase is irrelevant for a single isolated qubit. But when that qubit is part of a larger entangled system, what appears to be a “global” phase on the subsystem becomes a relative phase in the full system:
# Two-qubit state: (|00> + |11>) / sqrt(2) (Bell state)
bell = np.array([1, 0, 0, 1], dtype=complex) / np.sqrt(2)
# Apply a "global" phase of -1 to the second qubit only
# This means: |0> -> |0>, |1> -> -|1> on qubit 2
# Which is a Z gate on qubit 2
Z2 = np.kron(np.eye(2), np.array([[1, 0], [0, -1]]))
after_Z = Z2 @ bell
print(f"Before Z2: {bell}")
print(f"After Z2: {after_Z}")
# Before: [0.707, 0, 0, 0.707]
# After: [0.707, 0, 0, -0.707] -- different state! Relative phase changed.
# The probabilities in the computational basis are unchanged
print(f"Probs before: {np.abs(bell)**2}")
print(f"Probs after: {np.abs(after_Z)**2}")
# Both [0.5, 0, 0, 0.5], but entanglement correlations differ
Summary
Complex numbers are not a convenience in quantum computing. They are a necessity. The interference effects that give quantum algorithms their power require amplitudes that can cancel and reinforce, and that requires complex arithmetic. Real-valued amplitudes can represent some quantum states, but they cannot capture the full range of quantum dynamics.
The key points to remember:
- Amplitudes are complex numbers; probabilities are their squared magnitudes
- Complex multiplication is rotation and scaling, which is exactly what phase gates perform
- The full group of single-qubit unitaries (SU(2)) requires complex numbers; real matrices (SO(2)) are too restrictive
- Phase is the angle of a complex amplitude in the complex plane
- Global phase is unobservable for isolated qubits; relative phase drives interference
- Inner products require the complex conjugate on the left vector
- The Bloch sphere parameterization uses two real numbers from the complex amplitude ratio
- Phase kickback transfers eigenvalue phases to control qubits, enabling quantum phase estimation
- The quantum Fourier transform is built from roots of unity, complex numbers on the unit circle
- Density matrix off-diagonal entries encode coherence; decoherence destroys them
- Quantum state tomography recovers complex amplitudes through measurements in multiple bases
Next Steps
- Linear Algebra for Quantum Computing: vectors, matrices, and the full mathematical framework
- Bra-Ket Notation Explained: the notation that physicists use for complex quantum states
- How Quantum Algorithms Work, showing how interference is put to work in algorithm design
Was this tutorial helpful?