Decoding Quantum Errors on the Steane code with Belief Propagation & Catalyst

Symplectic Representation

To formalize stabilizers and errors, we use symplectic vectors. For \(n\) physical qubits, any \(n\)-qubit Pauli operator can be represented as a binary vector of length \(2n\):

where:

• \(v_i = 1\) if the Pauli has an X on qubit \(i\) (0 otherwise),

• \(u_i = 1\) if the Pauli has a Z on qubit \(i\) (0 otherwise),

• and \(v_i = u_i = 1\) if the Pauli has a Y on qubit \(i\), since \(X \cdot Z \propto Y\).

For example, the operator \(XZIY\) on 4 qubits corresponds to \((1,0,0, 1 | 0,1,0, 1)\). The stabilizer group is then generated by a set of such symplectic vectors, forming a stabilizer matrix of size \(m \times 2n\) (with \(m\) generators).

The commutation condition between two Pauli operators \((v|u)\) and \((v'|u')\) is captured by their symplectic inner product:

For a valid stabilizer code, all generators must commute, which implies that the symplectic inner product between any two rows of the stabilizer matrix must be zero.

Measurement and Syndrome

When an error \(e\) occurs, represented as a symplectic vector \((v_e | u_e)\), the syndrome is calculated by taking the symplectic inner product of \(e\) with each stabilizer generator. This produces a syndrome bit \(s_i\) for each generator:

where \((v_g^{(i)} | u_g^{(i)})\) is the \(i\)-th generator.

Below, we show the basic quantum circuit for extracting a syndrome value. The ancilla qubit is initialized in the \(|+\rangle\) state, and controlled operations are applied based on the stabilizer generators. Finally, the ancilla is measured in the \(X\)-basis to obtain the syndrome value. Check out Arthur Pesah’s excellent blog post series 1 on Stabilizer codes for a deeper introduction.

from typing import Callable, Optional, Dict, Union, Sequence

import pennylane as qml

syndromes = ["XXZIZIX", "XXIIZZI"]
dev = qml.device("lightning.qubit", wires := max(map(len, syndromes)) + 1)


@qml.qnode(device=dev)
def ancilla_assisted_syndrome_extraction(syndromes: list[str]):
    ancilla = wires - 1
    for i, syndrome in enumerate(syndromes):
        qml.Hadamard(ancilla)
        for i, s in enumerate(syndrome):
            if s == "X":
                qml.CNOT(wires=[ancilla, i])
            elif s == "Z":
                qml.CZ(wires=[ancilla, i])
        qml.Hadamard(ancilla)
        qml.measure(ancilla)
        qml.Barrier()
        ancilla += 1


print(qml.draw(ancilla_assisted_syndrome_extraction, show_all_wires=True)(syndromes))

0: ────╭X──────────────────────||────╭X───────────────────||─┤
1: ────│──╭X───────────────────||────│──╭X────────────────||─┤
2: ────│──│──╭Z────────────────||────│──│─────────────────||─┤
3: ────│──│──│─────────────────||────│──│─────────────────||─┤
4: ────│──│──│──╭Z─────────────||────│──│──╭Z─────────────||─┤
5: ────│──│──│──│──────────────||────│──│──│──╭Z──────────||─┤
6: ────│──│──│──│──╭X──────────||────│──│──│──│───────────||─┤
7: ──H─╰●─╰●─╰●─╰●─╰●──H──┤↗├──||────│──│──│──│───────────||─┤
8: ────────────────────────────||──H─╰●─╰●─╰●─╰●──H──┤↗├──||─┤

Decoding: The Classical Half of QEC

Once you have syndrome bits from your stabilizer measurements, you need to figure out what error likely occurred–this is the job of the decoder. Formally, given the syndrome, you’re solving for the most probable error, usually called the maximum likelihood estimate (MLE) for the error.

However, exact MLE decoding depends on the precise information of your noise model and is generally computationally intractable (NP-Hard) because \(n\) one-bit syndrome measurements can take on \(2^n\) unique values. In practice, we rely on approximate methods tuned to assumptions about the noise model.

CSS Codes: Simplifying the Structure

CSS (Calderbank–Shor–Steane) codes are a special class of stabilizer code where the generators are split into X-type and Z-type operators. Their symplectic vectors look like this:

• X-type generator: \((v | 0)\) (only Xs)

• Z-type generator: \((0 | u)\) (only Zs)

This allows us to represent the stabilizers with two \(m \times n\) matrices:

• \(H_X\) for X-type generators

• \(H_Z\) for Z-type generators

The commutation condition to ensure that all generators are simultaneously observable is:

which ensures that all X and Z stabilizers commute pairwise. When measuring syndromes:

• X-type stabilizers detect Z errors via \(s_X = H_X e_Z^T \pmod{2}\).

• Z-type stabilizers detect X errors via \(s_Z = H_Z e_X^T \pmod{2}\).

This separation makes decoding modular, allowing you to handle X and Z errors independently. When we introduce the Steane code later, you’ll see these matrices explicitly and how they simplify syndrome calculation and decoding. See a similar diagram below for the CSS code cycle structure.

The Sum-Product Algorithm

The BP decoder is based on the sum-product algorithm, which computes marginal probabilities over the graph. Here’s the procedure:

1. Initialization

   Each variable node \(v\) sends an initial message to its neighboring checks that reflects the intrinsic belief about whether an error has occurred. This is the log-likelihood ratio (LLR) based on the physical error rate \(p\):

   \[L_0 = \log\frac{1 - p}{p}\]

   This expresses the prior belief: if \(p\) is small (e.g., 0.01), then \(L_0\) is positive, favoring no error; if \(p\) is close to 0.5, \(L_0\) is near zero (no strong prior). In general \(p\) is a parameter of the algorithm that can be tuned to your specific noise source.

2. Variable-to-Check Update

   Each variable node updates its message to a neighboring check \(c\) by combining its intrinsic belief with the incoming messages from other connected checks:

   \[m_{v \to c} = L_0 + \sum_{c' \in N(v) \setminus c} m_{c' \to v}\]

   Here:

   • \(m_{v \to c}\) is the message from variable \(v\) to check \(c\).

   • \(N(v)\) is the set of checks connected to variable \(v\).

   • \(m_{c' \to v}\) are messages received from neighboring checks other than \(c\).

3. Check-to-Variable Update

   Each check node updates its message to a neighboring variable \(v\) based on the syndrome bit \(s_c\) and the incoming messages from the other variables connected to it:

   \[m_{c \to v} = (-1)^{s_c} \; 2 \, \operatorname{arctanh} \biggl( \prod_{v' \in N(c) \setminus v} \tanh\frac{m_{v' \to c}}{2} \biggr)\]

   Here:

   • \(m_{c \to v}\) is the message from check \(c\) to variable \(v\).

   • \(s_c\) is the syndrome bit for check \(c\) (0 if the stabilizer is satisfied, 1 if violated).

   • \(N(c)\) is the set of variables connected to check \(c\).

   • The \(\tanh\) and \(\operatorname{arctanh}\) functions implement the sum-product rule for combining binary parity checks derived from classical probability theory.

   What’s going on? This formula indicates that if the product of incoming \(\tanh\) terms is close to +1 or -1, it means there is a strong belief about whether the parity is satisfied or violated. The \(\operatorname{arctanh}\) converts that back into an LLR-style message. The \((-1)^{s_c}\) factor flips the sign if the syndrome is 1, signaling that a parity error was detected.

4. Iteration

   Steps 2 and 3 are repeated for a fixed number of iterations (e.g., 10–20) or until the messages converge (i.e., stop changing significantly). Traditional theory and heuristics in error correction say to repeat \(BP\) roughly on the order of \(O(n)\).

5. Decision Rule

   After the iterations, each variable node computes its posterior LLR by summing its intrinsic belief and all incoming messages:

   \[L_v = L_0 + \sum_{c \in N(v)} m_{c \to v}\]

   The decoder then makes a hard decision:

   • If \(L_v < 0\), it guesses \(e_v = 1\) (error detected).

   • If \(L_v > 0\), it guesses \(e_v = 0\) (no error).

Why This Works

Belief propagation is exact on tree-like graphs, where no cycles exist. However, even on Tanner graphs, which are never tree-like, it provides a good approximation to the maximum-likelihood decoder by using only local, iterative computations. Nevertheless, its performance can degrade when the Tanner graph contains many short cycles—a common characteristic of many popular quantum codes, which can lead to poor convergence. In practice, further extensions like BP-OSD 3, BP-LSD 4 or Ambiguity Clustering 5 are used to fix these issues.

See the following summary article 6 as well as Chapter 5 in Bayesian Reasoning and Machine Learning 7 for a deeper dive into message passing algorithms on graphs.

BP in JAX

Below, we implement a BP decoder using JAX broken down into it’s core components.

Before we can pass messages, we need to establish the connectivity between nodes. The _build_graph function scans the parity‑check matrix once and records, for every variable node, which checks touch it and vice versa. We convert the neighbour lists to tuples so they become immutable, hashable static data. JAX can then embed their values as compile‑time constants in the XLA program and reliably reuse the compiled kernel multiple times. A cool thing about JAX/XLA is that when using simple static parameters like the ones below, the individual integers it contains are baked into the XLA program as compile‑time constants, so we can truly compile a high performance decoder for our specific parity check matrix.

def _build_graph(
    pcm: ArrayLike,
) -> tuple[tuple[tuple[int, ...], ...], tuple[tuple[int, ...], ...]]:
    """
    Pre‑compute variable‑node and check‑node neighbors.

    Returns
    -------
    var_neighbors : tuple[tuple[int, ...], ...]  # length = n
    check_neighbors : tuple[tuple[int, ...], ...]  # length = m
    """
    m, n = pcm.shape
    vars_, checks_ = [[] for _ in range(n)], [[] for _ in range(m)]

    for c in range(m):
        for v in range(n):
            if pcm[c, v]:
                vars_[v].append(c)
                checks_[c].append(v)

    return tuple(map(tuple, vars_)), tuple(map(tuple, checks_))

A nice way to visulaize this Tanner graph is using the networkx package. Below is an example on the Steane code.

import matplotlib.pyplot as plt
import networkx as nx

vars, checks = _build_graph(H_steane)
G = nx.Graph()
num_vars = len(vars)
num_checks = len(checks)

# build the nx graph object from our vars and checks
for v in range(num_vars):
    G.add_node(f"v{v}", bipartite=0)
for c in range(num_checks):
    G.add_node(f"c{c}", bipartite=1)
for c in range(num_checks):
    for v in checks[c]:
        G.add_edge(f"c{c}", f"v{v}")

pos = nx.bipartite_layout(G, nodes=[f"v{i}" for i in range(num_vars)])

plt.figure(figsize=(10, 7))
nx.draw(G, pos, with_labels=True, node_color="skyblue", node_size=500, font_weight="bold")
plt.title("Bipartite Graph for H_steane", fontsize=16)
plt.show()

The _c2v_update helper function performs one full sweep of check‑to‑variable updates (step 3 of the sum‑product algorithm). It takes the previous messages, the syndrome, the neighbor tables, and two scalars (L_int for the intrinsic log‑likelihood ratio and eps for numerical safety). It loops only over existing edges, multiplies the relevant \(\operatorname{tanh}\) terms, clips the product, applies \(\operatorname{arctanh}\), and writes the new message into the next matrix.

def _c2v_update(
    m_c2v_prev: ArrayLike,
    syndrome: ArrayLike,
    var_nei: tuple[tuple[int, ...], ...],
    check_nei: tuple[tuple[int, ...], ...],
    L_int: float,
    eps: float,
) -> ArrayLike:
    """
    Compute the next round of check‑to‑variable messages.
    """
    m, n = m_c2v_prev.shape
    m_c2v_next = jnp.zeros_like(m_c2v_prev)

    # Loop over checks (outer) then their vars (inner)
    for c in range(m):
        Vc = check_nei[c]
        if len(Vc) < 2:
            continue  # degree‑1 checks carry no new info

        for v in Vc:
            prod = 1.0
            # product over all v' ≠ v in this check
            for v_p in Vc:
                if v_p == v:
                    continue
                incoming = L_int
                for c_p in var_nei[v_p]:
                    if c_p != c:
                        incoming += m_c2v_prev[c_p, v_p]
                prod *= jnp.tanh(incoming / 2.0)

            prod = jnp.clip(prod, -1.0 + eps, 1.0 - eps)
            msg = ((-1) ** syndrome[c]) * 2.0 * jnp.arctanh(prod)
            m_c2v_next = m_c2v_next.at[c, v].set(msg)

    return m_c2v_next

Once the main loop finishes, we still need a hard decision. The function _posterior_llrs folds every final check‑to‑variable message for bit v into its intrinsic LLR, yielding the posterior belief for that bit. A negative value means “error likely,” a positive value means “no error.”

def _posterior_llrs(
    m_c2v_final: ArrayLike, var_nei: tuple[tuple[int, ...], ...], L_int: float
) -> ArrayLike:
    """
    Combine intrinsic LLR with all incoming messages.
    """
    n = m_c2v_final.shape[1]
    llr = jnp.full(n, L_int)
    for v in range(n):
        for c in var_nei[v]:
            llr = llr.at[v].add(m_c2v_final[c, v])
    return llr

build_bp_decoder serves as the main entry point for compiling our decoder. It takes the parity‑check matrix and channel error rate, builds the parity graph, pre‑computes the intrinsic LLR, and returns a JIT‑compiled function _decode.

Inside _decode, the following steps are executed:

1. All messages are zero-initialized.

2. _c2v_update is called inside a jax.lax.fori_loop for max_iter rounds.

3. Final messages are converted to posterior LLRs with _posterior_llrs.

4. A binary error vector is output by thresholding the LLRs at zero.

Because the whole _decode body is wrapped in @jax.jit, the first call compiles everything into an XLA kernel; subsequent calls run at full device speed.

def build_bp_decoder(
    parity_check_matrix: ArrayLike,
    error_rate: float,
    max_iter: int = 10,
    epsilon: float = 1e-9,
) -> Callable[[ArrayLike], ArrayLike]:
    """
    Return a JIT‑compiled BP decoder for the given code and channel.

    Parameters
    ----------
    parity_check_matrix : array‑like (m, n)
    error_rate : float              # BSC crossover probability p
    max_iter : int
    epsilon : float                 # numerical safety margin
    """
    pcm = jnp.asarray(parity_check_matrix, dtype=jnp.int32)
    m, n = pcm.shape
    L_int = jnp.log((1.0 - error_rate) / error_rate)

    var_nei, check_nei = _build_graph(pcm)

    @jax.jit
    def _decode(syndrome: ArrayLike) -> ArrayLike:
        syndrome = jnp.asarray(syndrome, dtype=jnp.int32)

        # Initialise all messages to zero
        m_c2v = jnp.zeros((m, n), dtype=jnp.float32)

        # BP loop
        def body(_, msgs):
            return _c2v_update(msgs, syndrome, var_nei, check_nei, L_int, epsilon)

        m_c2v = jax.lax.fori_loop(0, max_iter, body, m_c2v)

        # Hard decision from posterior LLRs
        llr = _posterior_llrs(m_c2v, var_nei, L_int)
        return (llr < 0).astype(jnp.int32)

    # optionally we can force our decoder to compile right away by calling it on a test input
    _decode(jnp.zeros(m, dtype=jnp.int32))

    return _decode

Let’s test the performance of the BP decoder on the Steane code compared to the LUT decoder.

bp_steane = build_bp_decoder(H_steane, error_rate=0.05, max_iter=7)

n_bits = H_steane.shape[0]
correct = 0
total_syndromes = 2**n_bits

table_data = []
headers = ["Syndrome", "BP Estimated Error", "LUT Exact Error", "Match"]

for i in range(total_syndromes):
    syndrome_binary_string = f"{i:0{n_bits}b}"
    s_array = jnp.array([int(x) for x in syndrome_binary_string])

    # Get error patterns from BP decoder and LUT
    bp_pattern = bp_steane(s_array)
    lut_pattern = lut_steane(s_array)
    match = jnp.all(bp_pattern == lut_pattern)

    bp_pattern_str = "".join(map(str, bp_pattern.tolist()))
    lut_pattern_str = "".join(map(str, lut_pattern.tolist()))

    table_data.append([syndrome_binary_string, bp_pattern_str, lut_pattern_str, str(match)])

    # Increment correct count if patterns match
    if match:
        correct += 1

print(tabulate(table_data, headers=headers))

# Calculate and print the BP accuracy
accuracy = (correct / total_syndromes) * 100 if total_syndromes > 0 else 0
print(f"\nBP Accuracy: {accuracy:.2f}%")

Syndrome    BP Estimated Error    LUT Exact Error  Match
----------  --------------------  -----------------  -------
       000               0000000            0000000  True
       001               1000000            1000000  True
       010               0100000            0100000  True
       011               0010000            0010000  True
       100               0001000            0001000  True
       101               0000100            0000100  True
       110               0000010            0000010  True
       111               0000001            0000001  True

BP Accuracy: 100.00%

Before diving into the code, let’s test our belief‑propagation (BP) decoder on a bigger example: the n‑bit repetition code. This code stores each logical bit by repeating it \(n\) times (e.g. \(0 \mapsto 00\ldots0\) and \(1 \mapsto 11\ldots1\)). Its parity‑check matrix consists of \(n-1\) rows, each enforcing that two neighbouring bits are equal. Below, we measure how often the BP decoder corrects random errors on a 50‑bit repetition code and compare its success rate to an optimal maximum‑likelihood (ML) decoder, which simply picks the lower‑weight error pattern consistent with the observed syndrome.

def rep_code(n: int) -> ArrayLike:
    """
    Build the (n − 1) × n parity‑check matrix H for the [n, 1] repetition code.

    Each row enforces equality between two neighboring bits:
        H[i] has 1s in positions i and i+1, zeros elsewhere.
    """
    # First row: parity check on bits 0 and 1 → [1, 1, 0, 0, …, 0]
    first_row = jnp.zeros(n, dtype=jnp.int8).at[jnp.array([0, 1])].set(1)
    rows = [first_row]

    # Remaining rows: slide the two‑bit “window” to the right
    for _ in range(n - 2):
        rows.append(jnp.roll(rows[-1], 1))  # shift previous row by 1 position

    return jnp.stack(rows)  # shape = (n‑1, n)


@jax.jit
def ml_rep_decoder(syndrome: ArrayLike) -> ArrayLike:
    """
    Minimum‑weight decoder for the repetition code.

    Parameters
    ----------
    syndrome : ArrayLike, shape (n‑1,)
        The syndrome s = H e (mod 2).

    Returns
    -------
    error : ArrayLike, shape (n,)
        A lowest‑weight error vector consistent with `syndrome`.
    """
    # Candidate 1: assume e[0] = 0, then recover the rest via cumulative XOR.
    #   e[k+1] = e[k] ⊕ s[k]  ⇒  e = [0, cumsum(s) mod 2]
    e0 = jnp.concatenate((jnp.array([0], dtype=jnp.int32), jnp.mod(jnp.cumsum(syndrome), 2)))

    # Candidate 2: flip every bit (equivalent to choosing e[0] = 1).
    e1 = (e0 + 1) & 1  # fast “add‑one then mod 2”

    # Compare Hamming weights.
    w0, w1 = jnp.sum(e0), jnp.sum(e1)

    # Return the lighter candidate (ties resolved in favour of e0).
    return jax.lax.cond(w0 <= w1, lambda _: e0, lambda _: e1, operand=None)

We run a short experiment on a \(50\) bit repetition code. We sample 10,000 random syndromes vectors and compute the accuracy of our BP decoder compared to our baseline ml_rep_decoder

H_rep = rep_code(n := 50)
bp_rep = build_bp_decoder(parity_check_matrix=H_rep, error_rate=0.1, max_iter=n)

# sample random syndromes
N = 10_000
key = jax.random.PRNGKey(0)
syndromes = jax.random.randint(key, shape=(N, n - 1), minval=0, maxval=2)

# use jax to map the decoder over the syndromes
# since our decoders are jit compiled jax functions they can be used with jax.vmap
success_rate = jnp.mean(
    jnp.all(jax.vmap(ml_rep_decoder)(syndromes) == jax.vmap(bp_rep)(syndromes), axis=1)
)

print(f"Decoding success rate: {success_rate * 100:.2f}%")

Decoding success rate: 85.73%