Level 1 · Quantum Computing: Complex Numbers and Vectors

The memory of a quantum computer is not a row of bits. It is a list of numbers, one for every bit string the machine could show you, and those numbers are complex: arrows in a plane, each with a length and a direction. A program multiplies that list by matrices; reading the result turns the squared lengths of the arrows into probabilities. On this path you build the whole machine yourself, as a state-vector simulator in plain Python, and then run the famous quantum algorithms on it. No physics is needed: everything is complex numbers, vectors and matrices.

This level is the toolbox. The first half is complex numbers in Python: arithmetic, length and mirror image, angles, Euler's formula, and the fact that multiplying by an arrow of length 1 turns things. The second half is vectors and matrices with complex entries, kept as plain Python lists: the inner product with its all-important conjugate, length, the matrix-vector and matrix products, the conjugate transpose, matrices that keep every length, the Kronecker product and eigenvectors. Nothing here is quantum yet. Next level these lists become qubits.

The paths meet here. The Signals, Systems & Control path uses the same turning arrows for its Fourier transform at Level 3, and the quantum Fourier transform near the end of this path is its close cousin. The Machine Learning path starts with vectors as lists, dot products and lengths; here the same ideas gain a conjugate. Weights that add up to 1 are probabilities, the language of the Statistics path, and the bit strings and search problems of the Algorithms path come back as what a quantum computer works on.

Complex numbers in Python

A complex number z = x + iy has a real part x and an imaginary part y, where i is a number with i^2 = -1. Draw it as an arrow from the origin to the point (x, y) of the complex plane: the real part reaches right, the imaginary part up. Python has complex numbers built in and writes i as 1j (engineers write j for it, and Python followed them).

z = 2 + 3j
w = 1 - 1j
print(z + w, z - w, z * w, z / w)
print(z.real, z.imag, type(z).__name__)
print(1j * 1j, complex(2, 3) == z)

(3+2j) (1+4j) (5+1j) (-0.5+2.5j)
2.0 3.0 complex
(-1+0j) True

Adding works part by part, which puts the arrows tip to tail. Multiplying uses i^2 = -1: (2+3i)(1-i) = 2 - 2i + 3i - 3i^2 = 5 + i. Python divides as well. complex(x, y) builds a number from two floats, and .real and .imag are always floats.

Two more tools are needed all the time. The modulus |z| = \sqrt{x^2 + y^2} is the arrow's length, abs(z). The conjugate \bar z = x - iy is the mirror image in the real axis, z.conjugate(). Ints and floats have .conjugate() too (it changes nothing), so code that conjugates works on lists that mix them.

import math, cmath
z = 3 - 4j
print(abs(z), z.conjugate(), (5).conjugate(), (2.5).conjugate())
print(cmath.sqrt(-9))
try:
    math.sqrt(-9)
except ValueError:
    print("math.sqrt refuses negative numbers")

5.0 (3+4j) 5 2.5
3j
math.sqrt refuses negative numbers

The module math works with real numbers only; its twin cmath has the same functions for complex numbers.

Common mistakes: writing j or 4 j: a bare j is a variable name (a NameError) and 4 j is a syntax error; write 1j and 4j. z.conjugate without brackets, which is the method itself, not the mirrored number. Using math where cmath is needed.

Practice: Think of a Complex Number

The weight of an arrow

The modulus squared |z|^2 = x^2 + y^2 is the square of the arrow's length. Call it the arrow's weight: it is never negative, and it is the number this whole path keeps coming back to. Python offers three ways to compute it, and one tempting way that is wrong:

z = 1 + 3j
print(abs(z) ** 2, z.real ** 2 + z.imag ** 2)
print(z * z.conjugate(), z * z)

10.000000000000002 10.0
(10+0j) (-8+6j)

abs(z) ** 2 takes a square root and squares it again, so the last digit can wobble; z.real ** 2 + z.imag ** 2 avoids that. The product with the conjugate is real because (x + iy)(x - iy) = x^2 - i^2y^2 = x^2 + y^2:

Python still calls it a complex number, (10+0j), with imaginary part 0; take .real if you need a float. But z * z is not a weight at all: it is another complex number, and for 1j it is -1.

A list of arrows has a total weight, and dividing each weight by the total gives its share. The shares of any list that is not all zeros add up to 1, since dividing every weight by the same total makes the total 1. Do not confuse the weight of the list with the weight of the sum of the arrows:

arrows = [1 + 1j, 2, -1j]
total = sum(z.real ** 2 + z.imag ** 2 for z in arrows)
print("total weight:", total)
print("share of the arrow 2:", round(abs(arrows[1]) ** 2 / total, 4))
print("weight of the sum:", abs(sum(arrows)) ** 2)

total weight: 7.0
share of the arrow 2: 0.5714
weight of the sum: 9.0

The weights are 2, 4 and 1, so the arrow 2 holds 4/7 of the total. The arrows themselves add up to 3 + 0i, of weight 9: the imaginary parts of 1+1j and -1j cancelled. Adding arrows can cancel or reinforce, adding weights never cancels. That difference is what quantum algorithms exploit.

Aside (optional): where are the arrows? Nobody ever sees the arrows of a real quantum computer. Reading it out gives one outcome, and repeating the run gives each outcome about as often as its share says. The arrows are the bookkeeping that predicts those frequencies. Keeping them exact, inside a chip cooled to near absolute zero or in atoms held still by lasers, is what makes the machines so hard to build.

Common mistakes: z * z or z ** 2 for the weight. Forgetting to square abs(z). Dividing by the number of arrows instead of the total weight.

Practice: How Heavy Is Each Arrow?

Polar form: length and angle

Instead of its two coordinates, an arrow can be given by its length r = |z| and its angle \theta, measured anticlockwise from the positive real axis, in radians (a full turn is 2\pi, half a turn \pi). That is its polar form:

cmath.phase(z) returns the angle, cmath.polar(z) both numbers, and cmath.rect(r, theta) goes back:

import cmath, math
for z in [2j, -3, -4 - 3j, 1 - 1j]:
    angle = cmath.phase(z)
    print(z, abs(z), round(angle, 4), round(math.degrees(angle), 1))
print(cmath.polar(1 - 1j))
print(cmath.rect(2, math.pi / 6))
print(cmath.phase(-1), cmath.phase(complex(-1, -0.0)))

2j 2.0 1.5708 90.0
-3 3 3.1416 180.0
(-4-3j) 5.0 -2.4981 -143.1
(1-1j) 1.4142135623730951 -0.7854 -45.0
(1.4142135623730951, -0.7853981633974483)
(1.7320508075688774+0.9999999999999999j)
3.141592653589793 -3.141592653589793

cmath.phase always answers in the range (-\pi, \pi]: positive in the upper half of the plane, negative in the lower half, and exactly \pi for the negative real axis. The last line shows a trap: Python keeps a sign on zero, and an imaginary part of -0.0, which a multiplication can produce, counts as just below the axis, so the phase flips from \pi to -\pi.

Many tasks want angles in another range, for example [0, 2\pi) like a dial. Adding or subtracting whole turns of 2\pi never changes the direction, so any range can be reached. For [0, 2\pi) one % does it, because Python's % with a positive right side is never negative. For (-\pi, \pi], shift by whole turns until the angle is inside, and mind which end is included:

import math
def wrap(a):                      # into (-pi, pi]
    while a > math.pi:
        a -= 2 * math.pi
    while a <= -math.pi:
        a += 2 * math.pi
    return a

for a in [4.0, -2.0, 7.5, -math.pi]:
    print(a, "->", round(a % (2 * math.pi), 4), "and", round(wrap(a), 4))

4.0 -> 4.0 and -2.2832
-2.0 -> 4.2832 and -2.0
7.5 -> 1.2168 and 1.2168
-3.141592653589793 -> 3.1416 and 3.1416

Common mistakes: degrees where radians are expected (math.radians and math.degrees convert). Computing the angle as math.atan(y / x): it divides by zero on the imaginary axis and gives 1+1j and -1-1j the same angle; use cmath.phase(z) or math.atan2(y, x). Trusting the angle of 0, which has no direction at all.

Practice: Distance and Direction on a Dial, Which Way Did the Weather Vane Turn?, One More Arrow to Cancel the Hum

Euler's formula

The arrow of length 1 at angle \theta has a famous short name:

This is Euler's formula. It means the polar form can be written z = r\,e^{i\theta}, and Python computes it as r * cmath.exp(1j * theta):

import cmath, math
theta = math.pi / 3
print(cmath.exp(1j * theta))
print(complex(math.cos(theta), math.sin(theta)))
print(abs(cmath.exp(2.5j)))
print(cmath.exp(1j * math.pi))
try:
    math.exp(1j)
except TypeError:
    print("math.exp refuses complex numbers")

(0.5000000000000001+0.8660254037844386j)
(0.5000000000000001+0.8660254037844386j)
0.9999999999999999
(-1+1.2246467991473532e-16j)
math.exp refuses complex numbers

Every e^{i\theta} has length 1, up to rounding. Half a turn, e^{i\pi}, is -1 in exact maths; in floats a tiny imaginary part of about 10^{-16} is left over (more on that below). The 1j matters: cmath.exp(2.5) is the real number e^{2.5} \approx 12.18, a stretch, not a turn.

Common mistakes: math.exp(1j * theta), which raises an error. Forgetting the 1j. Passing degrees.

Practice: Paces and Bearings on a Treasure Map

Multiplying turns arrows

Write two numbers in polar form and multiply:

The lengths multiply and the angles add. So multiplying by a number of length 1, e^{i\varphi}, keeps every length and turns by \varphi: anticlockwise for a positive \varphi, clockwise for a negative one. Multiplying by 1j is a quarter turn anticlockwise, by -1 half a turn, by -1j a quarter turn clockwise:

import cmath, math
p = (2 * cmath.exp(0.5j)) * (3 * cmath.exp(1.2j))
print(round(abs(p), 6), round(cmath.phase(p), 6))
print((4 + 1j) * 1j, (4 + 1j) * -1, (4 + 1j) * -1j)
turn = cmath.exp(1j * math.radians(30))
print(round(abs((4 + 1j) * turn), 6), round(abs(4 + 1j), 6))

6.0 1.7
(-1+4j) (-4-1j) (1-4j)
4.123106 4.123106

Lengths 2 and 3 give 6, angles 0.5 and 1.2 give 1.7. The 30-degree turn leaves the length of 4+1j alone.

Powers are repeated products, so the n-th power multiplies the angle by n and raises the length to the power n. This is De Moivre's formula:

import cmath, math
step = cmath.exp(0.3j)
print(round((step ** 5).real, 6), round((step ** 5).imag, 6))
print(round(math.cos(1.5), 6), round(math.sin(1.5), 6))
points = [1.5 * step ** k for k in range(12)]
plot([z.real for z in points], [z.imag for z in points], kind="scatter", title="1.5, turned 0.3 rad at a time")

0.070737 0.997495
0.070737 0.997495

Five turns of 0.3 radians make 1.5 radians, so the fifth power lands on \cos 1.5 + i\sin 1.5. Press Run to see the twelve points walk round a circle of radius 1.5, 0.3 radians apart.

Common mistakes: forgetting that a factor whose length is not 1 also stretches. Turning the wrong way: a negative angle turns clockwise.

Turning about another point

Multiplying by e^{i\varphi} turns everything about the origin. To turn about another point c, the pin, move the picture so that the pin sits at the origin, turn, and move it back:

A door hinged at 1 + i has its handle at 3 + i. Open it by 60 degrees anticlockwise:

import cmath, math
hinge = 1 + 1j
handle = 3 + 1j
turn = cmath.exp(1j * math.pi / 3)
opened = (handle - hinge) * turn + hinge
print(round(opened.real, 4), round(opened.imag, 4))
print("distance to the hinge:", round(abs(handle - hinge), 6), "->", round(abs(opened - hinge), 6))
wrong = handle * turn
print("turned about 0 instead:", round(wrong.real, 4), round(wrong.imag, 4))

2.0 2.7321
distance to the hinge: 2.0 -> 2.0
turned about 0 instead: 0.634 3.0981

Seen from the hinge, the handle is at handle - hinge = 2; turned by 60 degrees that becomes 2e^{i\pi/3} = 1 + 1.7321i, and moving back adds the hinge. The handle stays 2 units from the hinge, as a door handle should. Turning the handle itself about 0 would rip the door off its hinge.

Common mistakes: multiplying p instead of p - c. Adding c back before multiplying. Using cmath.exp(angle) without the 1j.

Practice: Turning a Sketch About a Pin

Floats are almost never equal

A quarter turn is e^{2\pi i/4} = i, and four quarter turns make a full one. In floats:

import cmath, math
w = cmath.exp(2j * math.pi / 4)
print(w)
print(w ** 4)
print(w ** 4 == 1, abs(w ** 4 - 1) < 1e-9)

(6.123233995736766e-17+1j)
(1-2.4492935982947064e-16j)
False True

\pi is not exactly a float, so cos of a quarter turn is 6 \times 10^{-17} instead of 0, and the errors grow a little with every operation. So never compare floats with ==. Compare their distance with a small tolerance: abs(a - b) < tol. For complex numbers abs(a - b) is the distance between two points of the plane, so the same test works. For lists, test every entry with all(...).

To print complex results readably, round the two parts separately (round does not accept complex numbers). Adding 0.0 turns a -0.0 into 0.0:

import cmath, math
def show(z, digits=4):
    return complex(round(z.real, digits) + 0.0, round(z.imag, digits) + 0.0)

print(show(cmath.exp(2j * math.pi / 4)), show(cmath.exp(1j * math.pi)), show(complex(-1e-17, -2.71828)))

1j (-1+0j) -2.7183j

Rounding is for printing only. A problem that compares answers with a tolerance wants the unrounded values.

Common mistakes: == on floats. round(z) on a complex number (a TypeError). Rounding a result that is compared with a tolerance.

Roots of unity

Split the unit circle into n equal steps. The points e^{2\pi i k/n} for k = 0, 1, \dots, n-1 all satisfy w^n = 1, because by De Moivre \left(e^{2\pi i k/n}\right)^n = e^{2\pi i k}, a whole number of turns. They are the n-th roots of unity:

import cmath, math
def show(z, digits=4):
    return complex(round(z.real, digits) + 0.0, round(z.imag, digits) + 0.0)

n = 6
roots = [cmath.exp(2j * math.pi * k / n) for k in range(n)]
print([show(w) for w in roots])
print(all(abs(w ** n - 1) < 1e-9 for w in roots))
print("sum:", show(sum(roots)))
plot([w.real for w in roots], [w.imag for w in roots], kind="scatter", title="sixth roots of unity")

[(1+0j), (0.5+0.866j), (-0.5+0.866j), (-1+0j), (-0.5-0.866j), (0.5-0.866j)]
True
sum: 0j

Press Run to see a regular hexagon. The roots add up to 0, for every n \ge 2, and the reason is a turn: multiplying every root by w = e^{2\pi i/n} moves each one a step along to the next, so the set, and its sum S, stay the same. Then wS = S with w \ne 1, which leaves only S = 0. Arrows spread evenly round a circle cancel completely; this cancellation drives the quantum Fourier transform near the end of this path.

Common mistakes: testing w ** n == 1. Forgetting the 2 * math.pi, or the 1j.

Practice: The Clock Hand That Skips Ticks, Every Number Whose Power Is z

Vectors, and one factor for all

A vector here is a Python list of complex numbers (ints and floats allowed). Adding two vectors and scaling one by a number work entry by entry; zip walks two lists side by side:

u = [1 + 1j, 2, 3j]
v = [0.5, 1j, 1]
print([a + b for a, b in zip(u, v)])
print([2j * a for a in u])

[(1.5+1j), (2+1j), (1+3j)]
[(-2+2j), 4j, (-6+0j)]

Multiplying every entry by the same number c of length 1, c = e^{i\varphi}, turns all the arrows by the same angle. No weight changes, so no share changes either:

import cmath
a = [1 + 1j, 2.0, 3j]
c = cmath.exp(0.7j)
b = [c * x for x in a]
print([round(abs(x) ** 2, 6) for x in a])
print([round(abs(x) ** 2, 6) for x in b])
print(all(abs(y - c * x) < 1e-9 for x, y in zip(a, b)))
twisted = [a[0], c * a[1], a[2]]
print([round(abs(x) ** 2, 6) for x in twisted])

[2.0, 4.0, 9.0]
[2.0, 4.0, 9.0]
True
[2.0, 4.0, 9.0]

b and a are equal up to one factor of modulus 1: there is a single c with |c| = 1 and b_k = c\,a_k for every k. twisted has the same weights too, but it turned only its middle entry, so no single factor turns a into it: equal weights are not enough. A factor of another length, like 2 or 2j, is not allowed either; it would multiply every weight by 4. Later in this path, two states that differ only by one such factor count as the same state, and the next level gives that factor its name.

Common mistakes: comparing only the weights. Checking each entry with its own factor instead of one shared factor. Dividing b[k] / a[k] where a[k] is 0.

Practice: Same Snapshot, Turned Camera

The inner product

The inner product of two vectors of the same length multiplies matching entries and adds them up, after conjugating the entries of the first vector:

u = [1 + 2j, 1]
v = [3, 1j]
print(sum(a.conjugate() * b for a, b in zip(u, v)))
print(sum(a.conjugate() * b for a, b in zip(v, u)))
print(sum(a.conjugate() * b for a, b in zip(u, u)), sum(a * b for a, b in zip(u, u)))

(3-5j)
(3+5j)
(6+0j) (-2+4j)

\langle u|v\rangle = (1 - 2i)\cdot 3 + 1 \cdot i = 3 - 5i. Swapping the vectors conjugates the answer, \langle v|u\rangle = \overline{\langle u|v\rangle}, so which one carries the conjugate matters. The reason for the conjugate is the third line: a vector with itself gives

the total weight, real and never negative, the squared length of u. Without the conjugate u with itself gives -2+4j, which is no length at all. For real vectors the conjugate changes nothing, and this is the dot product of the Machine Learning path.

Two vectors with \langle u|v\rangle = 0 are orthogonal, the complex version of perpendicular. [1, 1j] and [1j, 1] look alike, but \overline{1}\cdot i + \overline{i}\cdot 1 = i - i = 0: they are orthogonal.

Common mistakes: leaving out the conjugate. Conjugating the second vector instead of the first. a.conjugate without brackets.

Practice: How Much Do Two Arrow Lists Agree?

Length and unit vectors

The length (or norm) of a vector is the square root of its total weight:

Dividing every entry by the length gives a unit vector, of length 1, pointing the same way. That is normalising. Every state of the simulator you will build is a unit vector, so its shares are just its weights.

import math
v = [1, 2j, 2]
length = math.sqrt(sum(abs(x) ** 2 for x in v))
unit = [x / length for x in v]
print(length)
print(unit)
print(math.sqrt(sum(abs(x) ** 2 for x in unit)))
print(sum(x * x for x in v))

3.0
[0.3333333333333333, 0.6666666666666666j, 0.6666666666666666]
1.0
(1+0j)

The weights are 1, 4 and 4, so the length is 3. The last line is a warning: squaring the entries themselves gives 1 - 4 + 4 = 1, a believable but wrong "length 1". Only the vector of zeros cannot be normalised: its length is 0 and it has no direction.

Common mistakes: x ** 2 instead of abs(x) ** 2. Forgetting the square root. Dividing by a length of 0.

Practice: Shrink It to Length One, Along the Line and Across It

Matrices and the matrix-vector product

A matrix is a list of rows: m[i] is row i and m[i][j] is the entry in row i, column j. A matrix acts on a vector by combining each row with the vector, entry by entry:

m = [[1, 1j],
     [0, 2],
     [1, -1]]
v = [2, 1j]
print(len(m), "rows of", len(m[0]), "entries; m[0][1] =", m[0][1])
for i, row in enumerate(m):
    entry = sum(row[j] * v[j] for j in range(len(v)))
    print("row", i, row, "->", entry)

3 rows of 2 entries; m[0][1] = 1j
row 0 [1, 1j] -> (1+0j)
row 1 [0, 2] -> 2j
row 2 [1, -1] -> (2-1j)

Row 0 gives 1 \cdot 2 + i \cdot i = 1. The answer has one entry per row: three here, although v has two. There is no conjugate in this product. A square matrix gives back as many entries as it takes in, and square matrices will be the gates of your simulator. The matrix [[0, -1], [1, 0]], for example, sends (x, y) to (-y, x), a quarter turn.

Common mistakes: walking down the columns instead of along the rows. Conjugating anything. Assuming the matrix is square.

Practice: The Mixing Desk

The matrix product, and why order matters

Applying B to a vector and then A to the result is the same as applying one matrix, the product AB:

Entry (i, j) combines row i of A with column j of B. The matrix applied first stands on the right, as in A(Bv) = (AB)v. Here A is a quarter turn and B doubles the second entry:

A = [[0, -1], [1, 0]]
B = [[1, 0], [0, 2]]
AB = [[sum(A[i][t] * B[t][j] for t in range(2)) for j in range(2)] for i in range(2)]
BA = [[sum(B[i][t] * A[t][j] for t in range(2)) for j in range(2)] for i in range(2)]
print("AB =", AB, " BA =", BA)
v = [3, 1]
Bv = [sum(B[i][j] * v[j] for j in range(2)) for i in range(2)]
print("B, then A:", [sum(A[i][j] * Bv[j] for j in range(2)) for i in range(2)])
print("AB at once:", [sum(AB[i][j] * v[j] for j in range(2)) for i in range(2)])

AB = [[0, -2], [1, 0]]  BA = [[0, -1], [2, 0]]
B, then A: [-2, 3]
AB at once: [-2, 3]

Stretch-then-turn and turn-then-stretch are different matrices: AB \ne BA in general. Like socks and shoes, order matters, and it will matter for gates too: the gates called X and Z next level have XZ \ne ZX. For matrices that are not square, the product of an r \times k and a k \times c matrix is r \times c: its rows come from A, its columns from B, and the sum over t runs over the k in the middle.

Common mistakes: multiplying entry by entry (A[i][j] * B[i][j]). Combining row i of A with row j of B. Writing the matrices in the order they are applied.

Practice: Socks Before Shoes

The dagger and unitary matrices

The conjugate transpose M^\dagger ("M dagger") swaps rows and columns and conjugates every entry:

zip(*m) hands out the columns of m, which become the rows of the dagger:

m = [[1, 2j, 0],
     [1 - 1j, 3, 4j]]
print(list(zip(*m)))
print([[x.conjugate() for x in col] for col in zip(*m)])

[(1, (1-1j)), (2j, 3), (0, 4j)]
[[1, (1+1j)], [-2j, 3], [0, -4j]]

Two rows of three become three rows of two. A square matrix U is unitary when

where I is the identity matrix, 1 on the diagonal and 0 elsewhere. Then U^\dagger undoes U, and U keeps the length of every vector: \|Uv\| = \|v\|. That is why every quantum gate is unitary. Here is one, checked with a tolerance, next to a matrix that is not unitary:

import math
s = 1 / math.sqrt(2)
U = [[s, -1j * s], [-1j * s, s]]
Ud = [[x.conjugate() for x in col] for col in zip(*U)]
UdU = [[sum(Ud[i][t] * U[t][j] for t in range(2)) for j in range(2)] for i in range(2)]
print(UdU)
print(all(abs(UdU[i][j] - (1 if i == j else 0)) < 1e-9 for i in range(2) for j in range(2)))

def length(v):
    return round(math.sqrt(sum(abs(x) ** 2 for x in v)), 6)
v = [1, 2]
for name, M in (("U", U), ("shear", [[1, 1], [0, 1]])):
    Mv = [sum(M[i][j] * v[j] for j in range(2)) for i in range(2)]
    print(name, length(v), "->", length(Mv))

[[(0.9999999999999998+0j), 0j], [0j, (0.9999999999999998+0j)]]
True
U 2.236068 -> 2.236068
shear 2.236068 -> 3.605551

The diagonal comes out as 0.9999999999999998, so == would wrongly say no. The shear [[1, 1], [0, 1]] turns [1, 2] into [3, 2], longer: it is not unitary.

Look at entry (i, j) of U^\dagger U: it is \sum_t \overline{U_{ti}}\,U_{tj}, the inner product of column i with column j. So U^\dagger U = I says exactly: every column has length 1, and any two different columns are orthogonal. Unitary is a word for square matrices only; a tall matrix can have perfect columns and still not be one.

Common mistakes: transposing without conjugating, or the other way round. Comparing with ==. Checking that the columns have length 1 but not that they are orthogonal.

Practice: Turn It Over, Then Reflect, Does the Machine Keep Every Length?

The Kronecker product

The Kronecker product (or tensor product) a \otimes b of a vector a of length m and a vector b of length n lists every product a[i] * b[j], mn of them, in a fixed order:

kron(a, b)[i * len(b) + j] = a[i] * b[j]

The index j runs fastest, like the last digit of a counter: first a[0] times all of b, then a[1] times all of b. That is the order of two nested loops with a outside and b inside, and also the order of a comprehension with two fors, where the first for is the outer loop:

print([s + t for s in "AB" for t in "xyz"])
a = [2, 1j]
b = [1, 3, 1 + 1j]
for i in range(len(a)):
    for j in range(len(b)):
        print(f"position {i * len(b) + j}: a[{i}] * b[{j}] = {a[i] * b[j]}")

['Ax', 'Ay', 'Az', 'Bx', 'By', 'Bz']
position 0: a[0] * b[0] = 2
position 1: a[0] * b[1] = 6
position 2: a[0] * b[2] = (2+2j)
position 3: a[1] * b[0] = 1j
position 4: a[1] * b[1] = 3j
position 5: a[1] * b[2] = (-1+1j)

There is no conjugate, and the order of the two vectors matters: b \otimes a holds the same six numbers in a different order. Weights multiply too, so the Kronecker product of two unit vectors is a unit vector. Your simulator will use it to join the states of separate parts into one state.

Matrices combine the same way. If B has p rows and q columns, kron(A, B)[i * p + k][j * q + l] = A[i][j] * B[k][l]. Pictured, A \otimes B is a grid of blocks, where block (i, j) is A[i][j] times all of B:

For A = \begin{pmatrix}1&2\\3&4\end{pmatrix} and B = \begin{pmatrix}0&1\\1&0\end{pmatrix} that is

and entry [2][1], with i = 1, k = 0, j = 0, l = 1, is A_{10}B_{01} = 3 \cdot 1 = 3.

Common mistakes: swapping the loops, which builds b \otimes a. Positions j * len(a) + i. A list of lists instead of one flat list. Conjugating.

Practice: Every Pair, in Order, A Matrix Made of Matrices

Eigenvectors: arrows a matrix only stretches

Most vectors change direction when a matrix acts on them. A nonzero vector v that the matrix only scales,

is an eigenvector of M, and the number \lambda (complex in general, possibly 0) is its eigenvalue. Take M = \begin{pmatrix}1 & i\\ 0 & 2\end{pmatrix} and test some vectors for \lambda = 2:

M = [[1, 1j], [0, 2]]
lam = 2
for v in ([1j, 1], [-1, 1j], [1, 1]):
    Mv = [sum(M[i][j] * v[j] for j in range(2)) for i in range(2)]
    ok = all(abs(Mv[i] - lam * v[i]) < 1e-9 for i in range(2))
    print(v, "->", Mv, ok)

[1j, 1] -> [2j, (2+0j)] True
[-1, 1j] -> [(-2+0j), 2j] True
[1, 1] -> [(1+1j), 2] False

[1j, 1] is doubled. So is [-1, 1j], which is 1j times [1j, 1]: any nonzero multiple of an eigenvector is an eigenvector with the same eigenvalue. Even after normalising to length 1 there is a free factor of modulus 1 left, so [0.7071j, 0.7071] and [-0.7071, 0.7071j] are the same answer, equal up to one factor of modulus 1. [1, 1] changes direction: it is not an eigenvector. The zero vector never counts, because M\,0 = \lambda\,0 for every \lambda.

To find the eigenvalues of a 2 by 2 matrix \begin{pmatrix}a&b\\c&d\end{pmatrix}, note that (M - \lambda I)v = 0 has a nonzero solution exactly when (a - \lambda)(d - \lambda) - bc = 0, a quadratic equation for \lambda. Here (1 - \lambda)(2 - \lambda) = 0 gives 1 and 2. Eigenvectors of gates will matter a great deal later in this path.

Common mistakes: reading \lambda off an entry where v is 0. Comparing with ==. Treating an eigenvalue of 0 as "no eigenvector" (if lam: is false for 0).

Practice: Arrows the Matrix Only Stretches, Two Directions That Stay Put

From formula to Python

Every formula of this level is a sum, and a sum is a loop, usually a generator expression inside sum. The index under \sum becomes the for, and the term becomes the expression in front of it. When two lists are walked side by side, zip names their entries; when a result is a whole list or matrix, a list comprehension (or two nested ones) holds one sum per entry.

• The weight |z|^2 = x^2 + y^2 is z.real ** 2 + z.imag ** 2, or abs(z) ** 2.
• The inner product \langle u|v\rangle = \sum_i \overline{u_i}\,v_i is sum(a.conjugate() * b for a, b in zip(u, v)).
• The length \|v\| = \sqrt{\sum_i |v_i|^2} is math.sqrt(sum(abs(x) ** 2 for x in v)).
• One entry of a matrix-vector product, (Mv)_i = \sum_j M_{ij}v_j, is sum(m[i][j] * v[j] for j in range(len(v))).
• One entry of a matrix product, (AB)_{ij} = \sum_t A_{it}B_{tj}, is sum(a[i][t] * b[t][j] for t in range(len(b))).
• One entry of the dagger, (M^\dagger)_{ij} = \overline{M_{ji}}, is m[j][i].conjugate(): note the swapped indices.
• The Kronecker rule (a \otimes b)_{i n + j} = a_i b_j is two loops, a outside, b inside.
• The sum of the roots of unity \sum_{k=0}^{n-1} e^{2\pi i k/n} is sum(cmath.exp(2j * math.pi * k / n) for k in range(n)).

All of them together check a rule that ties the dagger to the inner product, \langle u|Mv\rangle = \langle M^\dagger u|v\rangle, on some made-up numbers:

M = [[1, 2j], [1 + 1j, 3]]
u = [1j, 2]
v = [3, 1 - 1j]
Mv = [sum(M[i][j] * v[j] for j in range(2)) for i in range(2)]
Md = [[M[j][i].conjugate() for j in range(2)] for i in range(2)]
Mdu = [sum(Md[i][j] * u[j] for j in range(2)) for i in range(2)]
left = sum(a.conjugate() * b for a, b in zip(u, Mv))
right = sum(a.conjugate() * b for a, b in zip(Mdu, v))
print(left, right, abs(left - right) < 1e-9)

(14-5j) (14-5j) True

Reading the other way works too: when a description gives a formula, name each symbol in code first (M_{ij} is m[i][j], \overline{u_i} is u[i].conjugate()), then write the sum as a generator.

Next level these lists become qubits: a unit vector of two complex numbers is one qubit, a 2 by 2 unitary matrix is a gate, and the weights are the odds of what you see when you look.

Practice

These problems are required on this path at this level; the path adds one problem from each of the level's topics, picked for you and kept.

• How Heavy Is Each Arrow?
• Turning a Sketch About a Pin
• How Much Do Two Arrow Lists Agree?
• The Mixing Desk
• Does the Machine Keep Every Length?
• Every Pair, in Order

Checklist

You can now:

• Compute with Python's complex numbers (1j, .real, .imag, abs, .conjugate(), cmath) and avoid math where it refuses them.
• Compute the weight |z|^2 = z\bar z of an arrow, and the shares of a list of arrows, which add up to 1.
• Move between x + iy and polar form r e^{i\theta} with cmath.phase, cmath.polar, cmath.rect and cmath.exp, and bring an angle into [0, 2\pi) or (-\pi, \pi].
• Turn arrows by multiplying with e^{i\varphi}, about the origin or about any other point, and use De Moivre's formula for powers.
• Compare floats and complex numbers with a tolerance, never with ==, and print them rounded.
• List the n-th roots of unity and explain why they add up to 0.
• Tell when two lists are equal up to one factor of modulus 1.
• Compute inner products with the conjugate on the first vector, lengths, and unit vectors.
• Multiply matrices with vectors and with each other in the right order, take the dagger, and test whether a matrix is unitary.
• Build Kronecker products of vectors and matrices with the index rule, and check whether a vector is an eigenvector.