Get a correction ($e$) that satisfies $H.e=s$ from pymatching

Question

I am running some experiments on the surface code where I need the correction that satisfies $H.e=s$ and not the prediction of whether the observable has flipped or not from pymatching. I am trying to use the in-built API pymatching.Matching.from_check_matrix but this does not work when I have more than one round (rounds>1) in the circuit. Here is a code snippet of what I am doing:

circuit = stim.Circuit.generated('surface_code:rotated_memory_z', 
                                            after_clifford_depolarization=p,
                                            after_reset_flip_probability=p,
                                            before_measure_flip_probability=p,
                                            before_round_data_depolarization=p,
                                            distance=d,
                                            rounds=d)
dem = circuit.detector_error_model(decompose_errors=True)
sampler = dem.compile_sampler()
check_matrix, observable_matrix, priors = dem_to_check(dem)
matching = pymatching.Matching.from_check_matrix(check_matrix, error_probabilities=priors)

When using this code, Pymatching returns with the following error:

ValueError: `check_matrix` must contain at most two ones per column, but column 5 has 3 ones.

I am a little confused by this error, since the check matrix has been constructed from the DEM generated from a valid surface code circuit. What am I doing wrong here and how can I fix it?

score 2 · Accepted Answer · answered Mar 01 '25 at 15:16

I did some digging through ldpc=2.2.8, and I think I found your problem.

The main difference between the experiment with rounds=1 and rounds>1 is that some lines of the form error(0.01) D0 D1 ^ D2 appear in the latter cases while none appears in the former case. This means something has gone wrong in the handling of the decomposed errors.

We will focus on the following small snippet that fails in the exact same manner as yours:

import stim
import pymatching
from ldpc.ckt_noise import detector_error_model_to_check_matrices as dem_to_check
dem = stim.DetectorErrorModel("""
error(0.01) D0 D1 ^ D5 D6
error(0.02) D0 D2 ^ D5 D7
""")

dem2 = dem_to_check(dem)
check_matrix, observable_matrix, priors = dem2.check_matrix, dem2.observables_matrix, dem2.prior
matching = pymatching.Matching.from_check_matrix(check_matrix, error_probabilities=priors)

The return value of detector_error_model_to_check_matrices is a DemMatrices object from which you only extracted three fields: check_matrix, observable_matrix and priors.

By looking at how they are defined, we can learn that they handle the case where decomposed errors are not decomposed but remain as hyperedges (hence the "more than 2 ones per column" error from pymatching).

Two others fields exist and should give you the decomposed errors as edges: edge_check_matrix and edge_observables_matrix. However, merely using these fields is not enough:

import stim
import pymatching
from ldpc.ckt_noise import detector_error_model_to_check_matrices as dem_to_check
dem = stim.DetectorErrorModel("""
error(0.01) D0 D1 ^ D5 D6
error(0.02) D0 D2 ^ D5 D7
""")

dem2 = dem_to_check(dem)
check_matrix, observable_matrix, priors = dem2.edge_check_matrix, dem2.edge_observables_matrix, dem2.prior
matching = pymatching.Matching.from_check_matrix(check_matrix, error_probabilities=priors)

This now fails because the priors array is not coherent with the new matrices. There is no "edge_priors" in the DemMatrices class, so we will need to compute it ourselves using yet another of its field: hyperedge_to_edge_matrix.

import stim
import pymatching
from ldpc.ckt_noise import detector_error_model_to_check_matrices as dem_to_check
dem = stim.DetectorErrorModel("""
error(0.01) D0 D1 ^ D5 D6
error(0.02) D0 D2 ^ D5 D7
""")
dem2 = dem_to_check(dem)
check_matrix, observable_matrix, priors = dem2.edge_check_matrix, dem2.edge_observables_matrix, dem2.priors
edge_priors = dem2.hyperedge_to_edge_matrix @ priors
matching = pymatching.Matching.from_check_matrix(check_matrix, error_probabilities=edge_priors)

You can look at any sparse matrix with the todense method. From what I gathered, hyperedge_to_edge_matrix stores in its columns which edges correspond to each hyperedge, so I think the matrix product I used for edge_priors works to determine the correct prior probability of each edge.

Coming back to your stim experiment, it does not raise any error anymore, although I did not try to run it beyond this point:

import stim
import pymatching
from ldpc.ckt_noise import detector_error_model_to_check_matrices as dem_to_check
p = 0.01
d=3
circuit = stim.Circuit.generated('surface_code:rotated_memory_z', 
                                            after_clifford_depolarization=p,
                                            after_reset_flip_probability=p,
                                            before_measure_flip_probability=p,
                                            before_round_data_depolarization=p,
                                            distance=d,
                                            rounds=2)
dem = circuit.detector_error_model(decompose_errors=True)
dem2 = dem_to_check(dem)
check_matrix, observable_matrix, priors = dem2.edge_check_matrix, dem2.edge_observables_matrix, dem2.priors
edge_priors = dem2.hyperedge_to_edge_matrix @ priors
matching = pymatching.Matching.from_check_matrix(check_matrix, error_probabilities=edge_priors)

Get a correction ($e$) that satisfies $H.e=s$ from pymatching

1 Answers1