Level 6 · Quantum Computing: Noise and Error Correction

Every state so far was perfect: a unit vector that the gates turned exactly as asked. Real chips are not like that. A qubit leaks a little of itself into its surroundings, a gate turns slightly too far, a stray field flips a bit now and then. This level describes what that does and what can be done about it. The first half brings a new way to write a state, the density matrix, which can hold "this state or that one, with these odds" as well as a single state. With it come purity, the Bloch ball, the partial trace (what one part of a register looks like on its own), noise as channels with Kraus operators, and the fidelity that measures the damage.

The second half fights back. A bit can be protected by repeating it and taking a vote; a qubit cannot be copied, but it can be spread over three qubits so that one flipped wire is found and repaired without ever looking at the stored amplitudes. You will run the whole three-qubit code, see exactly how it fails, and turn it sideways to catch phase flips. The level ends with the variational eigensolver: a short circuit with a few dials, an "energy" made of Pauli readings, and a search for the dial setting that makes it smallest, the kind of algorithm today's noisy chips can already run.

The paths meet here. A mixture is a probability-weighted average, and the partial trace is a marginal of a joint table, both from the Statistics path; the majority vote fails with a binomial probability, and the noisy store is a seeded Monte Carlo simulation. A density matrix is Hermitian with non-negative eigenvalues, like the covariance matrix whose eigenvectors the Machine Learning path finds for PCA. Parity bits and XOR come from the Algorithms path, and its Level 6 black-box optimisation with a budget is exactly how the last problem is scored. numpy is available, and the examples use it for small full matrices.

A state written as a table

Level 2 showed that a global phase is invisible: [1j, 0] and [1, 0] are the same state. A cleaner way to write a state, with no leftover phase, is its density matrix, the outer product of the state with its own conjugate:

The conjugate sits on the column index j. In numpy that is np.outer(psi, psi.conj()). The examples of this level print matrices with a small helper show, which rounds and drops imaginary parts that are all zero:

import numpy as np
def show(a, d=4):
    a = np.round(np.asarray(a, dtype=complex), d) + 0
    return (a.real + 0).tolist() if not a.imag.any() else a.tolist()

psi = np.array([0.8, -0.6j])
rho = np.outer(psi, psi.conj())
print(show(rho))
print("rho[0][1] by hand:", 0.8 * (-0.6j).conjugate() + 0)
print("trace", round(float(np.trace(rho).real), 4), " Hermitian", np.allclose(rho, rho.conj().T))
turned = np.exp(0.7j) * psi                      # the same state with a global phase
print("same table:", np.allclose(np.outer(turned, turned.conj()), rho))
print("conjugate on the row instead:", show(np.outer(psi.conj(), psi)))

[[(0.64+0j), 0.48j], [-0.48j, (0.36+0j)]]
rho[0][1] by hand: 0.48j
trace 1.0  Hermitian True
same table: True
conjugate on the row instead: [[(0.64+0j), -0.48j], [0.48j, (0.36+0j)]]

Read the table like this. The diagonal holds the weights |\psi_i|^2, the Born-rule probabilities, so the trace (the sum of the diagonal, np.trace) is 1. The entries off the diagonal, the coherences, hold the relative phases. Swapping i and j conjugates an entry, \rho_{ji} = \overline{\rho_{ij}}, so \rho equals its own dagger: it is Hermitian. Trace 1 and Hermitian are two free checks on any density matrix you build. A global phase e^{i\phi} multiplies \psi_i by e^{i\phi} and \overline{\psi_j} by e^{-i\phi}, so it cancels in every entry.

Common mistakes: conjugating the row amplitude, which gives the transpose (the last line). Forgetting the conjugate, so \rho_{11} of [0, 1j] comes out as -1. np.outer(psi, psi) on a complex state.

Practice: The State Written as a Table

Mixtures: weigh the tables, never the amplitudes

A source that sends |\psi_k\rangle with probability p_k (a coin decides, and nobody tells you the coin) is a mixture. Its density matrix is the weighted sum of the tables:

Why this is enough: the chance of any outcome is the average, over the coin, of that outcome's chance for each |\psi_k\rangle, and each of those is a linear function of |\psi_k\rangle\langle\psi_k| (a diagonal entry, or a sum of entries after a gate). Averages of linear functions are linear functions of the average. So everything a measurement can tell you depends only on \rho, and two sources with the same \rho cannot be told apart by any experiment.

import numpy as np
def show(a, d=4):
    a = np.round(np.asarray(a, dtype=complex), d) + 0
    return (a.real + 0).tolist() if not a.imag.any() else a.tolist()
def dm(psi):
    psi = np.asarray(psi, dtype=complex)
    return np.outer(psi, psi.conj())

r = 1 / np.sqrt(2)
coin_z = 0.5 * dm([1, 0]) + 0.5 * dm([0, 1])              # |0> or |1>
coin_y = 0.5 * dm([r, 1j * r]) + 0.5 * dm([r, -1j * r])   # |+i> or |-i>
print(show(coin_y), np.allclose(coin_y, coin_z))
wrong = 0.5 * np.array([r, 1j * r]) + 0.5 * np.array([r, -1j * r])
print("mixing amplitudes instead:", show(wrong))
lopsided = 0.7 * dm([1, 0]) + 0.3 * dm([r, -r])           # |0> or |->
print(show(lopsided), round(float(np.trace(lopsided).real), 4))

[[0.5, 0.0], [0.0, 0.5]] True
mixing amplitudes instead: [0.7071, 0.0]
[[0.85, -0.15], [-0.15, 0.15]] 1.0

A fair coin between |0\rangle and |1\rangle and a fair coin between |{+i}\rangle and |{-i}\rangle give the same matrix, I/2 (and so does a fair coin between |+\rangle and |-\rangle): no experiment can tell those three sources apart. I/2 is the maximally mixed qubit, a perfect coin in every basis. Adding the amplitudes instead makes a superposition, a different thing: here it is not even normalised, and after normalising it is the pure state |0\rangle.

A matrix that is a single |\psi\rangle\langle\psi| is a pure state; one that is not is mixed. In code, build the sum in a fresh table (np.zeros((d, d), dtype=complex), or [[0j] * d for _ in range(d)] with plain lists).

Common mistakes: \sum_k p_k|\psi_k\rangle instead of \sum_k p_k|\psi_k\rangle\langle\psi_k|. [[0j] * d] * d, which repeats one row object so every row changes together. Weights that do not add up to 1.

Practice: Two Machines, One Output

Purity: how mixed is it?

One number says how far a matrix is from a single state, the purity

where \rho^2 is the matrix product \rho\rho and d is the size of the matrix. A pure state has \gamma = 1, because |\psi\rangle\langle\psi|\psi\rangle\langle\psi| = |\psi\rangle\langle\psi| (the middle \langle\psi|\psi\rangle is 1), so \rho^2 = \rho and its trace is 1. The other end is I/d, with \gamma = d \cdot \frac{1}{d^2} = \frac1d. Here a source sends |0\rangle or |+\rangle, with weight w on |+\rangle:

import numpy as np
def dm(psi):
    psi = np.asarray(psi, dtype=complex)
    return np.outer(psi, psi.conj())
def purity(rho):
    return float(np.trace(rho @ rho).real)

r = 1 / np.sqrt(2)
ws = [k / 20 for k in range(21)]
gs = [purity((1 - w) * dm([1, 0]) + w * dm([r, r])) for w in ws]
print([round(g, 4) for g in gs[::5]])
print(purity(np.eye(2) / 2), purity(np.eye(4) / 4), purity(dm([r, r])))
plot(ws, gs, title="purity of (1 - w)|0><0| + w|+><+|", xlabel="w", ylabel="Tr(rho^2)")

[1.0, 0.8125, 0.75, 0.8125, 1.0]
0.5 0.25 0.9999999999999996

Press Run to see the curve: 1 at both ends, where the source sends one state only, and lowest at w = \frac12. It stops at 0.75, not at \frac12: |0\rangle and |+\rangle overlap, so a coin between them is less uncertain than a coin between two orthogonal states. Note the last number: a pure state built from 1/\sqrt2 gives 0.9999999999999996, so "pure" must be tested with a tolerance, never with == 1.

Common mistakes: \sum_i \rho_{ii}^2, the diagonal squared, which forgets the coherences (it gives \frac12 for the pure |+\rangle). Multiplying entry by entry (rho * rho) and then taking the trace. Returning (1+0j) where a real number is wanted.

Practice: Is the Batch Pure?

What a measurement sees

Every reading comes out of \rho by a trace. Measuring in the computational basis gives outcome i with probability \rho_{ii}. The average reading of an observable P (Level 3's expectation value) is

For a pure state this is the old \langle\psi|P|\psi\rangle, since \mathrm{Tr}(P|\psi\rangle\langle\psi|) = \langle\psi|P|\psi\rangle; for a mixture it is the weighted average of the pure answers. For one qubit, the three Pauli expectations form the Bloch vector. Multiplying out the three traces once gives formulas in the entries:

A source that sends |{+i}\rangle or |0\rangle with even odds:

import numpy as np
def dm(psi):
    psi = np.asarray(psi, dtype=complex)
    return np.outer(psi, psi.conj())
X = np.array([[0, 1], [1, 0]])
Y = np.array([[0, -1j], [1j, 0]])
Z = np.diag([1, -1])

r = 1 / np.sqrt(2)
rho = 0.5 * dm([r, 1j * r]) + 0.5 * dm([1, 0])
by_trace = [round(float(np.trace(P @ rho).real), 4) + 0.0 for P in (X, Y, Z)]
by_entry = [round(float(v), 4) + 0.0 for v in (2 * rho[0, 1].real, -2 * rho[0, 1].imag, (rho[0, 0] - rho[1, 1]).real)]
print(by_trace, by_entry)
print("length", round(float(np.linalg.norm(by_trace)), 4))
print("P(0) =", round(float(rho[0, 0].real), 4), "= (1 + z) / 2")

[0.0, 0.5, 0.5] [0.0, 0.5, 0.5]
length 0.7071
P(0) = 0.75 = (1 + z) / 2

The Bloch vector of the mixture is the average of the two pure ones, (0, 1, 0) and (0, 0, 1), and the average of two points on a sphere lies inside it. So pure qubits sit on the surface of the Bloch sphere (length 1), mixed qubits inside the Bloch ball (length below 1), and I/2 at the centre. Each coordinate is the expectation of a \pm1 reading, which is why P(0) = (1 + z)/2. Going back, \rho = \frac12(I + xX + yY + zZ).

Common mistakes: the sign of y: |{+i}\rangle must give y = +1. np.trace(P * rho), entry by entry. Returning complex numbers with a zero imaginary part: take .real.

Practice: A Point Inside the Ball

Forgetting a qubit: the partial trace

Two stations share two qubits in the state \rho (4\times4). The first station holds qubit 0 only. What does it hold? Its own 2\times2 density matrix, the reduced density matrix, found by the partial trace over the qubit it cannot see.

Index i of two qubits is 2a + b, with a the bit of qubit 0 (the leftmost, most significant) and b that of qubit 1. Cut the 4\times4 matrix into four 2\times2 blocks: block (a, a') holds the rows with qubit 0 equal to a and the columns with qubit 0 equal to a'. Then:

• keep qubit 0: (\rho_0)_{aa'} = \sum_b \rho_{2a+b,\;2a'+b}, the trace of block (a, a');
• keep qubit 1: (\rho_1)_{bb'} = \sum_a \rho_{2a+b,\;2a+b'}, the sum of the two diagonal blocks.

In words: add up the entries where the forgotten qubit has the same value in the row and in the column. The diagonal of \rho_0 is then the marginal probability of qubit 0 (Statistics: add up a row of the joint table), and the off-diagonal terms survive only where the forgotten qubit agrees, because different values of a qubit nobody looks at can never interfere again. In general, when part A is the first qubits and part B the rest, every index splits as i = a\,d_B + b with d_B = 2^{\text{qubits in B}}, and the same rule holds. numpy's reshape does the splitting: rho.reshape(2, 2, 2, 2)[a, b, a2, b2] is rho[2 * a + b][2 * a2 + b2], and np.trace over a pair of axes adds the entries where those two indices agree.

The product of qubit 0 in 0.8|0\rangle + 0.6|1\rangle and qubit 1 in |{+i}\rangle should come apart into exactly those two:

import numpy as np
def show(a, d=4):
    a = np.round(np.asarray(a, dtype=complex), d) + 0
    return (a.real + 0).tolist() if not a.imag.any() else a.tolist()
def dm(psi):
    psi = np.asarray(psi, dtype=complex)
    return np.outer(psi, psi.conj())

r = 1 / np.sqrt(2)
rho = dm(np.kron([0.8, 0.6], [r, 1j * r]))       # qubit 0 first, as always
R = rho.reshape(2, 2, 2, 2)                       # R[a, b, a2, b2] = rho[2a + b][2a2 + b2]
keep0 = np.trace(R, axis1=1, axis2=3)             # b == b2: forget qubit 1
keep1 = np.trace(R, axis1=0, axis2=2)             # a == a2: forget qubit 0
print("keep qubit 0:", show(keep0))
print("by hand, entry [0][1]:", round(float((rho[0, 2] + rho[1, 3]).real), 4))
print("keep qubit 1:", show(keep1))

keep qubit 0: [[0.64, 0.48], [0.48, 0.36]]
by hand, entry [0][1]: 0.48
keep qubit 1: [[(0.5+0j), -0.5j], [0.5j, (0.5+0j)]]

Entry [0][1] of \rho_0 adds \rho_{0,2} (|00\rangle\langle10|) and \rho_{1,3} (|01\rangle\langle11|): qubit 1 is 0 in both row and column of the first, and 1 in both of the second. A corner such as \rho_{0,3}, where qubit 1 differs between row and column, never counts. Each answer has trace 1, a good check.

Common mistakes: swapping the two cases, so "keep qubit 0" returns qubit 1. Adding all four entries of a block. Treating qubit 0 as the rightmost bit, as some toolkits do.

Practice: Only One Half of the Pair, Keep These Qubits, Forget the Rest

Entanglement looks like noise

Now reduce an entangled pair, the Bell state |\Psi^-\rangle = (|01\rangle - |10\rangle)/\sqrt2, and compare it with a coin that picks |01\rangle or |10\rangle:

import numpy as np
def show(a, d=4):
    a = np.round(np.asarray(a, dtype=complex), d) + 0
    return (a.real + 0).tolist() if not a.imag.any() else a.tolist()
def dm(psi):
    psi = np.asarray(psi, dtype=complex)
    return np.outer(psi, psi.conj())
def halves(rho):
    R = rho.reshape(2, 2, 2, 2)
    return np.trace(R, axis1=1, axis2=3), np.trace(R, axis1=0, axis2=2)
def purity(rho):
    return round(float(np.trace(rho @ rho).real), 4)

r = 1 / np.sqrt(2)
bell = dm([0, r, -r, 0])
coin = 0.5 * dm([0, 1, 0, 0]) + 0.5 * dm([0, 0, 1, 0])
for name, rho in (("Bell pair", bell), ("coin", coin)):
    a, b = halves(rho)
    print(f"{name:9}  whole: purity {purity(rho)}   qubit 0: {show(a)}  qubit 1: {show(b)}")

leak = dm([0.6, 0, 0, 0.8])                       # 0.6|0> + 0.8|1> after a CNOT onto a qubit in |0>
print("before the leak:", show(dm([0.6, 0.8])), " after, qubit 0 alone:", show(halves(leak)[0]))

Bell pair  whole: purity 1.0   qubit 0: [[0.5, 0.0], [0.0, 0.5]]  qubit 1: [[0.5, 0.0], [0.0, 0.5]]
coin       whole: purity 0.5   qubit 0: [[0.5, 0.0], [0.0, 0.5]]  qubit 1: [[0.5, 0.0], [0.0, 0.5]]
before the leak: [[0.36, 0.48], [0.48, 0.64]]  after, qubit 0 alone: [[0.36, 0.0], [0.0, 0.64]]

The Bell pair is pure, yet each half on its own is I/2, the most mixed qubit there is. Locally it looks exactly like the coin: only a joint measurement (Level 3's correlations) tells them apart. That is also where noise comes from. In the last line a qubit 0.6|0\rangle + 0.8|1\rangle touched a second qubit, through a CNOT, and got entangled with it. Seen alone it keeps its probabilities, 0.36 and 0.64, but has lost its coherences: it is now a mixture. A real qubit does this, a little at a time, with the atoms and wires around it, and nobody keeps track of them.

Common mistakes: expecting the halves of a pure state to be pure. Concluding "no entanglement" because each half looks like noise; it is the other way round.

Eigenvalues: the mixture inside

A density matrix is Hermitian, so it has real eigenvalues and orthonormal eigenvectors (Level 1's eigenvectors), and np.linalg.eigvalsh finds the eigenvalues, smallest first. For a density matrix they are never negative and add up to 1: \rho is a mixture of its eigenvectors, with its eigenvalues as the odds. Two numbers follow. The purity is \sum_k \lambda_k^2, and the von Neumann entropy

counts the mixedness in bits: 0 for a pure state, 1 for I/2, n for the maximally mixed state of n qubits.

import math
import numpy as np
def entropy(rho):
    lam = np.linalg.eigvalsh(rho)
    return -sum(x * math.log2(x) for x in lam if x > 1e-12) + 0.0

rho = np.array([[0.75, -0.25j], [0.25j, 0.25]])     # the |+i>-or-|0> source from above
lam = np.linalg.eigvalsh(rho)
print(np.round(lam, 4).tolist(), round(float(sum(lam ** 2)), 4), round(float(np.trace(rho @ rho).real), 4))
print(round(entropy(rho), 4), entropy(np.eye(2) / 2), entropy(np.eye(8) / 8), entropy(np.diag([1.0, 0.0])))
print(entropy(np.outer([0.28, 0.96], [0.28, 0.96])))

[0.1464, 0.8536] 0.75 0.75
0.6009 1.0 3.0 0.0
0.0

The source sends |{+i}\rangle or |0\rangle, but the same \rho is also a mixture of two orthogonal states with odds 0.8536 and 0.1464; the purity agrees both ways. Eigenvalues that should be 0 can come out as -1e-17, which math.log2 rejects, so skip the tiny ones (the rule 0\log0 = 0). The same eigendecomposition \rho = V\,\mathrm{diag}(\lambda)\,V^\dagger (from np.linalg.eigh) gives functions of a matrix, such as its square root V\,\mathrm{diag}(\sqrt{\lambda})\,V^\dagger, which the fidelity of two mixed states needs.

Common mistakes: np.linalg.eig, which may return complex eigenvalues with tiny imaginary parts. Taking log2 of 0 or of rounding dust. Square-rooting the entries (np.sqrt(rho)) instead of the eigenvalues.

Practice: What the Two Halves Share

Channels: noise as Kraus operators

A gate is a unitary U, and on a density matrix it acts as \rho \mapsto U\rho\,U^\dagger. Noise is more general. A wire that applies X with probability q and leaves the qubit alone otherwise gives the mixture (1-q)\rho + q\,X\rho X^\dagger. Write the weights into the matrices, K_0 = \sqrt{1-q}\,I and K_1 = \sqrt q\,X, and it becomes

the general form of a noisy process, a channel, with its Kraus operators K_k. The square roots are there because each K appears twice. Not every list of matrices is a channel: the trace must stay 1. Since \mathrm{Tr}(K\rho K^\dagger) = \mathrm{Tr}(K^\dagger K\rho) (a trace does not change when the last factor moves to the front), that holds for every \rho exactly when

A unitary gate is the channel with one Kraus operator, since U^\dagger U = I. Here is the flip channel with q = 0.1, and the completeness check on a channel whose operators are not multiples of unitaries:

import numpy as np
def show(a, d=4):
    a = np.round(np.asarray(a, dtype=complex), d) + 0
    return (a.real + 0).tolist() if not a.imag.any() else a.tolist()
def dag(K):
    return K.conj().T

X = np.array([[0, 1], [1, 0]])
q = 0.1
flip = [np.sqrt(1 - q) * np.eye(2), np.sqrt(q) * X]
psi = np.array([0.8, 0.6])
rho = np.outer(psi, psi.conj())
out = sum(K @ rho @ dag(K) for K in flip)
print(show(rho), "->", show(out), round(float(np.trace(out).real), 4))
print("flip:   sum of K^dag K =", show(sum(dag(K) @ K for K in flip)))

g = 0.19
odd = [np.array([[1, 0], [0, np.sqrt(1 - g)]]), np.array([[0, np.sqrt(g)], [0, 0]])]
print("other:  sum of K^dag K =", show(sum(dag(K) @ K for K in odd)))
print("        sum of K K^dag =", show(sum(K @ dag(K) for K in odd)))

[[0.64, 0.48], [0.48, 0.36]] -> [[0.612, 0.48], [0.48, 0.388]] 1.0
flip:   sum of K^dag K = [[1.0, 0.0], [0.0, 1.0]]
other:  sum of K^dag K = [[1.0, 0.0], [0.0, 1.0]]
        sum of K K^dag = [[1.19, 0.0], [0.0, 0.81]]

The flip channel moved weight between the diagonal entries, 0.64 \to 0.9\cdot0.64 + 0.1\cdot0.36 = 0.612, and left the coherence 0.48 alone: for this real state X\rho X has the same off-diagonal entries as \rho. For the second channel the order matters: \sum K^\dagger K is the identity, so it is a valid channel, while \sum KK^\dagger is not. The test is always K^\dagger K, dagger on the left; the channel itself is K\rho K^\dagger, dagger on the right. Check with a tolerance, as for every float.

Common mistakes: K^\dagger\rho K. Testing \sum KK^\dagger, which rejects valid channels. A transpose without the conjugate, which breaks complex operators such as Y. Forgetting that \sqrt q, not q, goes into the operator.

Practice: The Wire That Smudges Qubits

Two standard noises

Two channels describe most of what happens to an idle qubit.

Depolarizing with probability p replaces the state by the maximally mixed one with probability p:

On a qubit it shrinks the Bloch vector by the factor 1 - p, straight toward the centre. It has a Kraus form too, \sqrt{1 - 3p/4}\,I and \sqrt{p/4} times each of X, Y, Z, but the formula above is simpler.

Amplitude damping with strength \gamma lets |1\rangle decay to |0\rangle:

K_1 moves |1\rangle to |0\rangle with weight \gamma; K_0 keeps |0\rangle and shrinks what is left of |1\rangle. That is the channel of the last section, with \gamma = 0.19.

import numpy as np
def show(a, d=4):
    a = np.round(np.asarray(a, dtype=complex), d) + 0
    return (a.real + 0).tolist() if not a.imag.any() else a.tolist()
def dm(psi):
    psi = np.asarray(psi, dtype=complex)
    return np.outer(psi, psi.conj())
def channel(rho, kraus):
    return sum(K @ rho @ K.conj().T for K in kraus)
def depolarize(rho, p):
    return (1 - p) * rho + p * np.eye(len(rho)) / len(rho)
def damping(g):
    return [np.array([[1, 0], [0, np.sqrt(1 - g)]]), np.array([[0, np.sqrt(g)], [0, 0]])]

p = 0.2
paulis = [np.eye(2), np.array([[0, 1], [1, 0]]), np.array([[0, -1j], [1j, 0]]), np.diag([1, -1])]
kraus = [np.sqrt(1 - 3 * p / 4) * paulis[0]] + [np.sqrt(p / 4) * P for P in paulis[1:]]
rho = dm([0.8, -0.6j])
print(show(depolarize(rho, p)), np.allclose(channel(rho, kraus), depolarize(rho, p)))

r = 1 / np.sqrt(2)
for name, psi in (("|1>", [0, 1]), ("|0>", [1, 0]), ("|+>", [r, r])):
    print(name, show(channel(dm(psi), damping(0.19))))

[[(0.612+0j), 0.384j], [-0.384j, (0.388+0j)]] True
|1> [[0.19, 0.0], [0.0, 0.81]]
|0> [[1.0, 0.0], [0.0, 0.0]]
|+> [[0.595, 0.45], [0.45, 0.405]]

|1\rangle keeps 0.81 of its weight and gives 0.19 to |0\rangle; |0\rangle is untouched, so damping drives every state toward |0\rangle, the north pole, rather than toward the centre. On |+\rangle it also shrinks the coherence, by \sqrt{1-\gamma} = 0.9.

Aside (optional): T1 and T2. On hardware, amplitude damping is energy loss: the state |1\rangle of a superconducting qubit holds a little more energy than |0\rangle and leaks it into the chip. Its time scale is called T1. The loss of coherences alone, as in the leak through a CNOT above, is dephasing, with time scale T2. Both are measured in microseconds to milliseconds today, and a gate takes tens of nanoseconds, so a circuit has room for some thousands of gates before the qubits have forgotten what they held.

Common mistakes: adding p\,I instead of p\,I/d. Applying depolarizing to the state vector. Putting \gamma instead of \sqrt\gamma into K_1.

Noise, step after step

On a chip the noise acts again at every clock tick, and when two channels act per tick, a problem must say in which order. The order matters whenever the channels do not commute. One tick on |1\rangle, with damping \gamma = 0.2 and a flip channel q = 0.1:

import numpy as np
def channel(rho, kraus):
    return sum(K @ rho @ K.conj().T for K in kraus)
def damping(g):
    return [np.array([[1, 0], [0, np.sqrt(1 - g)]]), np.array([[0, np.sqrt(g)], [0, 0]])]
def flip(q):
    return [np.sqrt(1 - q) * np.eye(2), np.sqrt(q) * np.array([[0, 1], [1, 0]])]

one = np.diag([0.0, 1.0])
a = channel(channel(one, damping(0.2)), flip(0.1))
b = channel(channel(one, flip(0.1)), damping(0.2))
print("damping, then flip: P(1) =", round(float(a[1, 1].real), 4))
print("flip, then damping: P(1) =", round(float(b[1, 1].real), 4))

rho, p1, gam = one, [], []
for tick in range(41):
    p1.append(float(rho[1, 1].real))
    gam.append(float(np.trace(rho @ rho).real))
    rho = channel(channel(rho, damping(0.2)), flip(0.1))
print([round(v, 4) for v in p1[:6]], round(p1[-1], 4), round(gam[-1], 4))
plot(p1, label="P(1)", title="|1> under damping 0.2, then flips 0.1, per tick", xlabel="tick")
plot(gam, label="purity")

damping, then flip: P(1) = 0.74
flip, then damping: P(1) = 0.72
[1.0, 0.74, 0.5736, 0.4671, 0.3989, 0.3553] 0.2778 0.5988

Different orders, different answers (0.74 against 0.72). Over many ticks the qubit forgets where it started: press Run to see P(1) fall from 1 and level off at 0.2778 = 0.1/0.36, where the weight that damping takes from |1\rangle each tick equals the weight that the flips put back, the same balance as a steady state in the Signals & Control path. The purity dips and recovers a little, because the resting state leans toward |0\rangle. Keep the current \rho in one variable, and record before the first tick if the start is wanted too.

Common mistakes: applying the channels in the other order. Recording only after each tick, which loses the start. Building rho from the original state again every tick.

Practice: The Qubit That Forgets

Fidelity: how much damage?

To score the damage, compare what arrived, \rho, with the pure state that was sent, |\psi\rangle. The fidelity is

the probability that a test "is it |\psi\rangle?" says yes. No square root is taken. F = 1 exactly when \rho = |\psi\rangle\langle\psi|, and for a pure \rho = |\varphi\rangle\langle\varphi| it is |\langle\psi|\varphi\rangle|^2. When \rho is |\psi\rangle after a channel, put \rho = \sum_k K_k|\psi\rangle\langle\psi|K_k^\dagger into the formula: each term is \langle\psi|K_k|\psi\rangle\,\langle\psi|K_k^\dagger|\psi\rangle, a number times its conjugate, so

and the matrix \rho is never needed. In a register, noise often hits one qubit only. The plain route builds the operator on qubit q of n as a Kronecker product, qubit 0 first (Level 3). Three qubits in 0.6|000\rangle + 0.8|111\rangle, with amplitude damping \gamma = 0.19 on qubit 2:

import numpy as np
def on_qubit(K, q, n):                               # K on qubit q of n, qubit 0 leftmost
    return np.kron(np.kron(np.eye(2 ** q), K), np.eye(2 ** (n - 1 - q)))

g = 0.19
damp = [np.array([[1, 0], [0, np.sqrt(1 - g)]]), np.array([[0, np.sqrt(g)], [0, 0]])]
psi = np.zeros(8, dtype=complex)
psi[0], psi[7] = 0.6, 0.8
big = [on_qubit(K, 2, 3) for K in damp]
rho = sum(K @ np.outer(psi, psi.conj()) @ K.conj().T for K in big)
print("from rho:      ", round(float(np.vdot(psi, rho @ psi).real), 6))
print("sum of overlaps:", round(sum(abs(np.vdot(psi, K @ psi)) ** 2 for K in big), 6))
plus = np.array([1, 1]) / np.sqrt(2)
ov = [np.vdot(plus, K @ plus) for K in damp]
print("|+> alone:", round(sum(abs(o) ** 2 for o in ov), 6), "  overlaps added first:", round(abs(sum(ov)) ** 2, 6))

from rho:       0.876096
sum of overlaps: 0.876096
|+> alone: 0.95   overlaps added first: 1.364095

K_0 keeps |000\rangle and shrinks |111\rangle to 0.9 \cdot 0.8, an overlap of 0.36 + 0.576 = 0.936; K_1 moves |111\rangle to |110\rangle, which has no overlap with |\psi\rangle. So F = 0.936^2 = 0.876096. For a lone |+\rangle under the same damping both operators overlap: the right answer is 0.95, while adding the overlaps before squaring gives more than 1, which no fidelity can be. For a register of many qubits, Level 3's pairs of amplitudes apply K to one qubit without the 2^n\times2^n matrix.

Common mistakes: taking a square root. Squaring the sum of overlaps instead of summing the squares. Putting K on the wrong end of the Kronecker product. Forgetting the conjugate on \psi (np.vdot takes care of it).

Practice: Damage in Transit, How Alike Are Two Blurred States?

Repetition and the majority vote

Before qubits, bits. A wire flips each bit with probability p, independently. Send the bit three times and let the receiver take the majority: the vote is wrong when two or three copies flip. Counting all eight flip patterns for p = 0.2:

from itertools import product
p = 0.2
fail = 0.0
for pattern in product([0, 1], repeat=3):            # 1 = this copy flipped
    prob = 1.0
    for f in pattern:
        prob *= p if f else 1 - p
    if sum(pattern) >= 2:
        fail += prob
print(round(fail, 6), round(3 * p ** 2 - 2 * p ** 3, 6))
ps = [k / 50 for k in range(51)]
plot(ps, ps, label="one copy", title="probability that the received bit is wrong", xlabel="p")
plot(ps, [3 * x ** 2 - 2 * x ** 3 for x in ps], label="majority of three")

0.104 0.104

Three patterns have two flips, each with probability p^2(1-p), and one has three: 3p^2(1-p) + p^3 = 3p^2 - 2p^3. For p = 0.2 that is 0.104, about half the error of a single copy; for p = 0.01 it is about 0.0003. Press Run: the curves cross at p = \frac12. Below it repetition helps, above it hurts, because then the wire is more often wrong than right and the vote trusts the wrong side. More copies help more, with binomial counts (Statistics path) in place of the hand count.

Common mistakes: counting one flip as a failure. Forgetting that two flips can happen in three ways. Expecting repetition to help when p > \frac12.

Practice: Outvoting a Noisy Wire

Spreading a qubit over three wires

A qubit cannot be copied (Level 3's aside on no-cloning), so a|0\rangle + b|1\rangle cannot become three copies of itself. It can be spread out: put it on qubit 0, two helper qubits in |0\rangle, and apply CNOT(0, 1) then CNOT(0, 2). Each CNOT copies the basis label of qubit 0, so |000\rangle stays and |100\rangle becomes |111\rangle:

That is the three-qubit bit-flip code. The damage it guards against is an X on one wire, which flips that qubit's bit in every basis index; on qubit k of n the bit has the mask 1 << (n - 1 - k), since qubit 0 is the most significant. With full matrices, for a = 0.8, b = -0.6:

import numpy as np
n = 3
def cnot(c, t):
    U = np.zeros((2 ** n, 2 ** n))
    for i in range(2 ** n):
        j = i ^ (1 << (n - 1 - t)) if (i >> (n - 1 - c)) & 1 else i
        U[j, i] = 1
    return U
def occupied(v):
    return [(format(i, "03b"), round(float(v[i].real), 4)) for i in range(len(v)) if abs(v[i]) > 1e-12]

a, b = 0.8, -0.6
start = np.kron([a, b], [1, 0, 0, 0])              # qubit 0 holds the qubit, qubits 1 and 2 start in |00>
code = cnot(0, 2) @ cnot(0, 1) @ start
print(occupied(code))
copies = np.kron(np.kron([a, b], [a, b]), [a, b])
print(int(np.sum(np.abs(copies) > 1e-12)), "nonzero amplitudes in three copies")
mask = 1 << (n - 1 - 1)                            # X on qubit 1
print(mask, occupied(np.array([code[i ^ mask] for i in range(8)])))

[('000', 0.8), ('111', -0.6)]
8 nonzero amplitudes in three copies
2 [('010', 0.8), ('101', -0.6)]

The code state has two nonzero amplitudes, three copies would have eight. An X on qubit 1 (mask 2) moves the amplitude of |000\rangle to |010\rangle and that of |111\rangle to |101\rangle: X never mixes amplitudes, it only moves them, so new[i] = old[i ^ mask]. A Z on qubit k instead negates every amplitude whose index has that bit set.

Common mistakes: taking the Kronecker product of three copies. 1 << k as the mask, which counts qubits from the right. Forgetting that the encoder runs on the register, not on the lone qubit.

Practice: One Qubit Spread Over Three Wires

Parity checks: the syndrome

How do you find the flipped wire without measuring the qubits, which would destroy a and b? Measure only parities. The check Z_0Z_1 asks "do qubits 0 and 1 agree?", the check Z_1Z_2 "do qubits 1 and 2 agree?". Their bits are

read from the bits b_0b_1b_2 of a basis state, and the pair (s_1, s_2) is the syndrome. The trick: in a|000\rangle + b|111\rangle both occupied basis states give (0, 0), and an X on one wire flips the same bit in both, so they still agree with each other. A measurement whose answer is the same for every occupied basis state is certain, and it leaves the state as it was: it learns where the damage is and nothing about a and b. In code, read the bits of any occupied index, the one with the largest amplitude being the safest. One row of the table, an X on qubit 2:

code = [0.8, 0, 0, 0, 0, 0, 0, -0.6]             # 0.8|000> - 0.6|111>
n, k = 3, 2
mask = 1 << (n - 1 - k)
damaged = [code[i ^ mask] for i in range(8)]
for i in range(8):
    if abs(damaged[i]) > 1e-9:
        b0, b1, b2 = (i >> 2) & 1, (i >> 1) & 1, i & 1
        print(format(i, "03b"), damaged[i], "syndrome", (b0 ^ b1, b1 ^ b2))

001 0.8 syndrome (0, 1)
110 -0.6 syndrome (0, 1)

A flip of qubit 2 changes only the second parity: syndrome (0, 1). Work out the other rows yourself, for no flip and for a flip of qubit 0 or qubit 1: the four syndromes come out different, so each one names its own repair. The repair applies X to the qubit that the syndrome blames; X undoes itself.

Common mistakes: reading the bits of index 0 when its amplitude is 0. Returning the parities as \pm1 or as booleans when bits are asked for. Measuring the qubits themselves, which collapses the stored state.

Practice: Two Parity Lamps on Three Wires, Find the Flipped Wire and Flip It Back

Decoding, and how the code fails

After the repair, run the encoder backwards: CNOT(0, 1) and CNOT(0, 2) again (each CNOT is its own inverse). The helpers return to |00\rangle and qubit 0 holds the result, at indices 0 (|000\rangle) and 4 (|100\rangle). Whatever the damage, the decoded qubit is the original |\psi\rangle with one of four things done to it, up to a global phase:

The code assumes at most one bit flip. With one, the syndrome is right and the result is |\psi\rangle. With two, the majority is wrong and the repair makes things worse; follow the amplitudes through once to see which row of the table you land on. Either way, the code fails exactly when the majority vote fails, so for independent flips with probability p its logical error rate is the 3p^2 - 2p^3 of the majority-vote section. And a Z on any wire gives a|000\rangle - b|111\rangle: both occupied states still have parities (0, 0), so the code sees nothing, repairs nothing, and the sign reaches the decoded qubit:

import numpy as np
n = 3
def cnot(c, t):
    U = np.zeros((2 ** n, 2 ** n))
    for i in range(2 ** n):
        j = i ^ (1 << (n - 1 - t)) if (i >> (n - 1 - c)) & 1 else i
        U[j, i] = 1
    return U
def z_on(v, k):
    m = 1 << (n - 1 - k)
    return np.array([-v[i] if i & m else v[i] for i in range(2 ** n)])
def same_state(u, v):                              # equal up to a global phase
    return bool(abs(abs(np.vdot(u, v)) - 1) < 1e-9)

psi = np.array([0.8, -0.6])
enc = cnot(0, 2) @ cnot(0, 1)
damaged = z_on(enc @ np.kron(psi, [1, 0, 0, 0]), 0)
i = int(np.argmax(np.abs(damaged)))
print("syndrome", (((i >> 2) ^ (i >> 1)) & 1, ((i >> 1) ^ i) & 1))
out = enc @ damaged                                 # decoding: the same two CNOTs
qubit = out[[0, 4]]
X, Z = np.array([[0, 1], [1, 0]]), np.diag([1, -1])
for name, M in (("I", np.eye(2)), ("X", X), ("Z", Z), ("XZ", X @ Z)):
    print(name, same_state(qubit, M @ psi))

syndrome (0, 0)
I False
X False
Z True
XZ False

One Z, and the stored qubit comes back as Z|\psi\rangle: for this code a phase flip is invisible. (Here enc doubles as the decoder because CNOT(0, 1) and CNOT(0, 2) commute and each undoes itself.) "Both" means X after Z, the matrix product XZ; since ZX = -XZ, the other order is the same state.

Common mistakes: reading qubit 0 at indices 0 and 1, which belong to qubit 2. Decoding before the repair. Calling Z|\psi\rangle "the same state": it differs by a relative sign, not a global one.

Practice: When Two Wires Flip at Once, A Thousand Nights in the Noisy Store

Turning the code sideways: phase flips

Level 2 showed that H swaps the roles of X and Z: HZH = X and HXH = Z. So a phase flip, seen through Hadamards, is a bit flip. The phase-flip code stores a|0\rangle + b|1\rangle as

the bit-flip code state with H on every qubit. A Z on one of its wires, followed by H on every qubit, becomes an X on the same wire of a|000\rangle + b|111\rangle, which the parity checks find. To repair: H on every qubit, syndrome and X as before, H on every qubit again.

import numpy as np
H = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
X, Z = np.array([[0, 1], [1, 0]]), np.diag([1, -1])
print(np.allclose(H @ Z @ H, X), np.allclose(H @ X @ H, Z))

H3 = np.kron(np.kron(H, H), H)
a, b = 0.8, -0.6
bitcode = np.zeros(8)
bitcode[0], bitcode[7] = a, b
phasecode = H3 @ bitcode                            # a|+++> + b|--->
Z1 = np.kron(np.kron(np.eye(2), Z), np.eye(2))      # Z on qubit 1
seen = H3 @ Z1 @ phasecode
print([(format(i, "03b"), round(float(v), 4)) for i, v in enumerate(seen) if abs(v) > 1e-12])

True True
[('010', 0.8), ('101', -0.6)]

After the H layer the Z on qubit 1 shows up as the bit-flip code state with qubit 1 flipped, 0.8|010\rangle - 0.6|101\rangle, and from there the bit-flip code's table takes over. This code is blind to X instead. Shor's nine-qubit code nests the two (three blocks of three) and catches either kind of error on any one qubit.

Aside (optional): why error correction is hard to build. Real codes, such as the surface code, spread one logical qubit over hundreds or thousands of physical ones, and measure their parity checks over and over while the computation runs, using extra qubits so that the data qubits are never measured directly. It only pays off below a threshold: if each physical part fails less often than roughly one time in a hundred to a thousand, a bigger code makes the logical qubit better; above it, a bigger code only adds more places to fail, like the vote above with p > \frac12. Chips first crossed that threshold in the mid-2020s.

Common mistakes: reading the syndrome from the raw phase-code state, where all eight basis states are occupied. Forgetting the final H layer. Expecting the phase-flip code to catch an X.

Practice: The Same Code, Turned Sideways

Pauli strings and their averages

A Pauli string such as "XYI" names one gate per qubit, the leftmost letter for qubit 0. On a basis state every letter does something simple: I nothing, X flips the bit, Z multiplies by -1 when the bit is 1, and Y flips the bit and multiplies by a factor, because Y = \begin{pmatrix}0&-i\\i&0\end{pmatrix}:

So a whole string sends each basis state to one other basis state, times \pm1 or \pm i: a signed permutation, never a mix. Its matrix (the Kronecker product of the letters, qubit 0 first) has a single nonzero entry in every column:

import numpy as np
L = {"I": np.eye(2), "X": np.array([[0, 1], [1, 0]]), "Y": np.array([[0, -1j], [1j, 0]]), "Z": np.diag([1, -1])}
def pauli(s):
    M = np.array([[1]])
    for ch in s:                                    # leftmost letter = qubit 0
        M = np.kron(M, L[ch])
    return M

P = pauli("XY")
for i in range(4):
    j = int(np.argmax(np.abs(P[:, i])))
    print(format(i, "02b"), "->", complex(P[j, i]), "*", format(j, "02b"))

psi = np.array([0.6, 0, 0, 0.8j])
terms = [(1.0, "ZZ"), (0.5, "XY"), (-0.3, "IX")]
vals = [round(float(np.vdot(psi, pauli(s) @ psi).real), 4) + 0.0 for c, s in terms]
print(vals, round(sum(c * v for (c, s), v in zip(terms, vals)), 4))

00 -> 1j * 11
01 -> -1j * 10
10 -> 1j * 01
11 -> -1j * 00
[1.0, 0.96, 0.0] 1.48

"XY" on |01\rangle: X flips qubit 0, and Y meets a 1 on qubit 1, so the result is -i\,|10\rangle. The expectation \langle\psi|P|\psi\rangle of a string is a real number between -1 and 1, the average of its \pm1 readings. A weighted sum of strings, H = \sum_k c_k P_k, is scored by

called the energy of the state by tradition, though no physics is needed: it is a weighted sum of averages. Here \langle ZZ\rangle = 1 (both occupied states have equal bits), \langle XY\rangle = 0.96 and \langle IX\rangle = 0, so the energy is 0.48 + 1 = 1.48. Building 2^n\times2^n matrices is fine for a few qubits; for more, use the signed permutation on each basis index directly.

Common mistakes: applying the letters from the right end of the index. The sign of Y: Y|0\rangle = +i|1\rangle. Forgetting the conjugate of \langle\psi|. Returning a complex number: take .real.

Practice: Adding Up the Pauli Readings

The variational idea: dials and an energy meter

Which state has the smallest energy? Its value is the smallest eigenvalue of the matrix H = \sum_k c_k P_k, which np.linalg.eigvalsh(M)[0] gives for a small matrix. For 50 qubits the matrix is far too big, but a chip can prepare a state and estimate its energy. The variational eigensolver (VQE) uses that: a short circuit with a few angles, the ansatz, prepares a state; the chip measures its energy; a classical program turns the dials and tries again, looking for the lowest reading.

The smallest ansatz has one qubit and one dial, |\psi(\theta)\rangle = R_y(\theta)|0\rangle = (\cos\frac\theta2, \sin\frac\theta2). Its three Pauli averages follow from the double-angle formulas:

the last because both amplitudes are real. The energy is a curve in \theta. A scan for H = 0.4\,I - 0.5\,Z + 1.2\,X, and one for H = Z + Y:

import math
import numpy as np
L = {"I": np.eye(2), "X": np.array([[0, 1], [1, 0]]), "Y": np.array([[0, -1j], [1j, 0]]), "Z": np.diag([1, -1])}
def energy(terms, theta):
    psi = np.array([math.cos(theta / 2), math.sin(theta / 2)])
    return sum(c * float(np.vdot(psi, L[s] @ psi).real) for c, s in terms)

ts = [k * 2 * math.pi / 120 for k in range(121)]
for terms in ([(0.4, "I"), (-0.5, "Z"), (1.2, "X")], [(1.0, "Z"), (1.0, "Y")]):
    es = [energy(terms, t) for t in ts]
    k = min(range(121), key=lambda i: es[i])
    M = sum(c * L[s] for c, s in terms)
    print(f"best on the grid: theta = {ts[k]:.4f}, energy {es[k]:.4f};  smallest eigenvalue {np.linalg.eigvalsh(M)[0]:.4f}")
    plot(ts, es, label=" + ".join(f"{c}{s}" for c, s in terms), title="energy of Ry(theta)|0>", xlabel="theta")

best on the grid: theta = 5.1313, energy -0.8996;  smallest eigenvalue -0.9000
best on the grid: theta = 3.1416, energy -1.0000;  smallest eigenvalue -1.4142

Press Run to see two waves. For the first Hamiltonian a grid of 121 angles comes within a thousandth of the smallest eigenvalue, -0.9; finding the exact bottom of such a curve, and the exact angle, is better done with a little trigonometry than with a finer grid. For Z + Y the best the dial can do is -1, but the smallest eigenvalue is -\sqrt2 \approx -1.4142: the state that reaches it has a complex amplitude, and R_y(\theta)|0\rangle never does. An ansatz can only find what it can reach. Choosing one that is short enough to run on a noisy chip, yet able to reach the answer, is the art of VQE.

On a chip each energy reading costs real time (thousands of shots), so a solver gets a budget of readings, and the problem with four dials counts them: one call too many raises an error. It is scored like the black-box problems of the Algorithms path's Level 6: to pass, your final energy must be at least as low as the best of budget random dial settings, and the quality measures, on a logarithmic scale, how much of the gap between that random search and the smallest eigenvalue you closed, every factor of ten counting the same. A meter that counts its calls is a few lines:

import math, random
def make_meter(f, budget):
    used = [0]
    def meter(theta):
        if used[0] >= budget:
            raise RuntimeError("budget spent")
        used[0] += 1
        return f(theta)
    return meter, used

f = lambda t: 0.4 - 0.5 * math.cos(t) + 1.2 * math.sin(t)
meter, used = make_meter(f, 20)
rng = random.Random(1)
best = min(meter(rng.uniform(0, 2 * math.pi)) for _ in range(20))
print(round(best, 4), used[0])
try:
    meter(0.0)
except RuntimeError as e:
    print("call 21:", e)

-0.8865 20
call 21: budget spent

Twenty random readings land within about 0.015 of the bottom here, because one dial is easy to cover; with four dials random guesses spread thin, and a search that uses what each reading says about the shape of the curve does far better. Stop by counting readings, never by the clock.

Common mistakes: using \theta/2 where \theta belongs (\langle Z\rangle = \cos\theta). An angle outside the asked range, such as -3\pi/4 instead of 5\pi/4. Expecting a real-amplitude ansatz to reach an eigenvalue that needs complex amplitudes. Spending one reading more than the budget.

Practice: One Dial, Lowest Reading, Four Dials and an Energy Meter

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.

• Only One Half of the Pair
• Is the Batch Pure?
• The Wire That Smudges Qubits
• Damage in Transit
• When Two Wires Flip at Once
• One Dial, Lowest Reading

Checklist

You can now:

• Build the density matrix \rho_{ij} = \psi_i\overline{\psi_j} of a state and of a mixture, check it (trace 1, Hermitian), and explain why a global phase disappears and why different mixtures can share one \rho.
• Compute the purity \mathrm{Tr}(\rho^2), between 1/d and 1, and tell pure from mixed with a tolerance.
• Read probabilities and expectations \mathrm{Tr}(P\rho) from \rho, and place a qubit in the Bloch ball with x = 2\,\mathrm{Re}\,\rho_{01}, y = -2\,\mathrm{Im}\,\rho_{01}, z = \rho_{00} - \rho_{11}.
• Take the partial trace of a two-qubit (or two-part) state with qubit 0 as the most significant bit, and explain why half of a Bell pair looks like I/2.
• Apply a channel \rho \mapsto \sum_k K_k\rho K_k^\dagger, check completeness with \sum_k K_k^\dagger K_k = I, and use depolarizing and amplitude-damping noise, in a stated order, over many steps.
• Measure damage with the fidelity F = \langle\psi|\rho|\psi\rangle = \sum_k|\langle\psi|K_k|\psi\rangle|^2, with the noise on one qubit of a register.
• Encode a qubit in the three-qubit bit-flip code, apply X and Z errors with the mask 1 << (n - 1 - k), read the syndrome from two parity bits, repair and decode.
• Explain why two flips and a single Z defeat the bit-flip code, name the four logical outcomes, and turn the code sideways with Hadamards to catch phase flips.
• Apply Pauli strings as signed permutations, with Y|0\rangle = i|1\rangle, and add up an energy \sum_k c_k\langle P_k\rangle.
• Describe the variational eigensolver: an ansatz with dials, an energy to minimise under a budget, and the smallest eigenvalue as the target it may or may not reach.