`predictions` and `fault_ids` in Stim and PyMatching for surface code decoding

Question

When I tried to use Stim and PyMatching to decode errors in a Rotated Planar code, I noticed that, Regardedless of the code distance, predictions.shape, as well as the shape of observable_flips, is always (n_shots, 1),

So I refered to the docs, which says the second dimension parameter is num_fault_ids.

So what is exactly the fault_id? matching.py offers such explanations:

    fault_ids: set[int] or int, optional
        The indices of any self-inverse faults which are flipped when the edge is flipped, and which should be tracked.
        This could correspond to the IDs of physical Pauli errors that occur when this
        edge flips (physical frame changes). Alternatively,
        this attribute can be used to store the IDs of any logical observables that are
        flipped when an error occurs on an edge (logical frame changes). In earlier versions of PyMatching, this
        attribute was instead named `qubit_id` (since for CSS codes and physical frame changes, there can be
        a one-to-one correspondence between each fault ID and physical qubit ID). For backward
        compatibility, `qubit_id` can still be used instead of `fault_ids` as a keyword argument.
        By default None

However these explanations brings me more questions. As I have learned, the decoder is fed with syndrome information (from stabilizer measurement results), and yields the most probable error pattern (which physical qubit has suffered from X/Z errors.) But the PyMatching decoder just gives a 0/1, and compares it with observable_flips.

So here is my question:

1. What is the meaning of fault_id in surface codes? And why num_fault_ids is always 1 in rotated codes? Is it an observable flip or something?
2. If it corresponds to an observable flip, then how does the PyMatching gets it? The decoding algorithm should only use the decoding graph to tell us which qubits went wrong, how does it predict whether there is a logical change? And why it can be used to evaluate the correctness of the decoder?

Here is an example with a simple distance-3 rotated planar code and only one shot.

import stim
import pymatching
import numpy as np  

def decode_error(circuit:stim.Circuit):
    n_shots = 1
    # sample circuit.
    sampler = circuit.compile_detector_sampler()
    detection_events, observable_flips = sampler.sample(n_shots, separate_observables=True)
    
    # configure decoder
    detector_error_model = circuit.detector_error_model(decompose_errors=True)
    matcher = pymatching.Matching.from_detector_error_model(detector_error_model)
    print(detector_error_model)
    print(matcher.edges())

    # run decoder
    predictions = matcher.decode_batch(detection_events)

    print(detection_events)
    print("predictions shape:",predictions.shape)
    print("predictions=",predictions)
    print(observable_flips)

    for shot in range(n_shots):
        if np.array_equal(predictions[shot], observable_flips[shot]):
            print("decode correct")

        
if __name__=='__main__':
    d=3
    round=d*d;  
    noise=0.1
    
    circuit = stim.Circuit.generated(
        "surface_code:rotated_memory_z",
        rounds=round,
        distance=d,
        before_round_data_depolarization=noise)

    decode_error(circuit)

and it gets such outputs:

predictions shape: (1, 1)
predictions= [[0]]

For a single round code, the matcher.edges() is:

[(0, None, {'fault_ids': set(), 'weight': 2.6390573296152584, 'error_probability': 0.06666666666666667}), (0, 1, {'fault_ids': set(), 'weight': 2.6390573296152584, 'error_probability': 0.06666666666666667}), (1, 2, {'fault_ids': set(), 'weight': 2.6390573296152584, 'error_probability': 0.06666666666666667}), (1, None, {'fault_ids': {0}, 'weight': 1.9509992185627845, 'error_probability': 0.12444444444444444}), (2, None, {'fault_ids': set(), 'weight': 1.9509992185627845, 'error_probability': 0.12444444444444444}), (2, 3, {'fault_ids': set(), 'weight': 2.6390573296152584, 'error_probability': 0.06666666666666667}), (3, None, {'fault_ids': {0}, 'weight': 2.6390573296152584, 'error_probability': 0.06666666666666667})]

Answer 1

In an X memory experiment, the X observable is either flipped or not flipped. There's only one logical qubit, regardless of the code distance. So there's always one bit of output. If there were two observables being evaluated, there'd be two bits of output. Output scales with observables, not distance.

If you want the individual prediction for the presence of every single edge of the matching graph, call pymatching.Matching.decode_to_edges_array instead of pymatching.Matching.decode_batch. Note that this will be slower (like maybe 3x slower? I forget exactly how much). pymatching has to a bunch of additional work to refine the prediction into individual edges.

Note that, in more complex decoders like tensor contraction decoders, there really isn't a sense in which an individual edge is present or not. Those kinds of decoders sum over enormous numbers of different possible edge configurations consistent with the detection events; not just the most likely one. It is a limitation of matching that it ultimately only considers one edge set, and you should be cautious of overly relying on this property by asking for that edge set.