pymatching for rotated surface code

Question

I am trying to do MWPM with pymatching for a rotated surface code with distance $L$, but I can't figure out how to initialize the Matching object in a simple way. I only want a single round of decoding. For $L=3$, I have 9 qubits with the following numbering scheme (X marks the checkerboard square acted on by the X stabilizers)

         0 --- 1 --- 2
         |  X  |     | X
         3 --- 4 --- 5
       X |     |  X  |       
         6 --- 7 --- 8

Following the tutorial provided for the toric code, I generate a CSS-style parity check matrix for these stabilizers:

[[1 1 0 1 1 0 0 0 0]
 [0 0 0 0 1 1 0 1 1]
 [0 0 1 0 0 1 0 0 0]
 [0 0 0 1 0 0 1 0 0]]

And then generate an example error and try to decode using pymatching:

matching = Matching(Hx)
noise = np.array([0, 0, 0, 1, 0, 0, 0, 0, 0]) # Z error on qubit '3'
s = Hx @ noise % 2 # correctly evaluates to [1, 0, 0, 1]
c = matching.decode(s) # predicts [0, 0, 0, 0, 0, 0, 0, 0, 0]

The prediction (regardless of what logical operators I use) doesn't see the error. As far as I can tell, pymatching is generating a graph that is mishandling the boundaries of the code. What did I do wrong, and what is the correct way to set up MWPM with pymatching for this rotated surface code?

Also I would appreciate any additional tutorial material for pymatching - without stim preferably :) I don't want/need to model circuit-level behavior or ancillary qubits.

AG47 · Accepted Answer · 2025-02-22T16:18:51.853

Using the following code, with the version 2.1.0 of pymatching and 1.25.2 of numpy, I get the expected result:

from pymatching import Matching
import numpy as np
import matplotlib.pyplot as plt
Hx = np.array([[1, 1, 0, 1, 1, 0, 0, 0, 0],
               [0, 0, 0, 0, 1, 1, 0, 1, 1],
               [0, 0, 1, 0, 0, 1, 0, 0, 0],
               [0, 0, 0, 1, 0, 0, 1, 0, 0]])
matching = Matching(Hx)
noise = np.array([0, 0, 0, 1, 0, 0, 0, 0, 0])
s = Hx @ noise % 2 # correctly evaluates to [1, 0, 0, 1]
c = matching.decode(s) # correctly predicts [0, 0, 0, 1, 0, 0, 0, 0, 0]
matching.draw()
plt.show()

There is nothing different from the code you posted, so I am not sure what exactly is wrong with it. As you mentionned in a comment, pymatching is not compatible with numpy>=2, so combining the two is a very likely explanation of the result discrepancy.

The matching graph used by pymatching is the following:

On this graph, nodes 0,1,2 and 3 are the four $X$ stabilizers and node 4 is the boundary node. Edges corresponds error mechanisms that flip neighboring nodes, so the non-equivalent $Z$ Pauli errors that could occur. There are only seven of them, because a $Z$ Pauli error on qubit 0 or qubit 1 (similarly qubit 7 or 8) has the same effect up to a $Z$-stabilizer of weight 2, so only one edge is used each time (edge 0 and 7, but no edge 1 and 8).

If you were to specify to pymatching flip probabilities, the weight of these edges would probably be higher than the other edges, because more error mechanisms correspond to them. Since you only give the parity-check matrix, all the edges are by default equiprobable (see matching.edges()).

pymatching for rotated surface code

1 Answers1