Data-driven decoding of quantum error correcting codes using graph
  neural networks
To leverage the full potential of quantum error-correcting stabilizer codes
it is crucial to have an efficient and accurate decoder. Accurate, maximum
likelihood, decoders are computationally very expensive whereas decoders based
on more efficient algorithms give sub-optimal performance. In addition, the
accuracy will depend on the quality of models and estimates of error rates for
idling qubits, gates, measurements, and resets, and will typically assume
symmetric error channels. In this work, instead, we explore a model-free,
data-driven, approach to decoding, using a graph neural network (GNN). The
decoding problem is formulated as a graph classification task in which a set of
stabilizer measurements is mapped to an annotated detector graph for which the
neural network predicts the most likely logical error class. We show that the
GNN-based decoder can outperform a matching decoder for circuit level noise on
the surface code given only simulated experimental data, even if the matching
decoder is given full information of the underlying error model. Although
training is computationally demanding, inference is fast and scales
approximately linearly with the space-time volume of the code. We also find
that we can use large, but more limited, datasets of real experimental data
[Google Quantum AI, Nature {\bf 614}, 676 (2023)] for the repetition code,
giving decoding accuracies that are on par with minimum weight perfect
matching. The results show that a purely data-driven approach to decoding may
be a viable future option for practical quantum error correction, which is
competitive in terms of speed, accuracy, and versatility.
Moritz Lange, Pontus Havström, Basudha Srivastava, Valdemar Bergentall, Karl Hammar, Olivia Heuts, Evert van Nieuwenburg, Mats Granath
# Data-driven decoding of quantum error correcting codes using graph neural networks

###### Abstract

To leverage the full potential of quantum error-correcting stabilizer codes it is crucial to have an efficient and accurate decoder. Accurate, maximum likelihood, decoders are computationally very expensive whereas decoders based on more efficient algorithms give sub-optimal performance. In addition, the accuracy will depend on the quality of models and estimates of error rates for idling qubits, gates, measurements, and resets, and will typically assume symmetric error channels. In this work, instead, we explore a model-free, data-driven, approach to decoding, using a graph neural network (GNN). The decoding problem is formulated as a graph classification task in which a set of stabilizer measurements is mapped to an annotated detector graph for which the neural network predicts the most likely logical error class. We show that the GNN-based decoder can outperform a matching decoder for circuit level noise on the surface code given only simulated experimental data, even if the matching decoder is given full information of the underlying error model. Although training is computationally demanding, inference is fast and scales approximately linearly with the space-time volume of the code. We also find that we can use large, but more limited, datasets of real experimental data [Google Quantum AI, Nature **614**, 676 (2023)] for the repetition code, giving decoding accuracies that are on par with minimum weight perfect matching. The results show that a purely data-driven approach to decoding may be a viable future option for practical quantum error correction, which is competitive in terms of speed, accuracy, and versatility.

## I Introduction

Quantum Error Correction (QEC) is foreseen to be a vital component in the development of practical quantum computing [1; 2; 3; 4; 5]. The need for QEC arises due to the susceptibility of quantum information to noise, which can rapidly accumulate and corrupt the final output. Unlike noise mitigation schemes where errors are reduced by classical post-processing [6; 7; 8], QEC methods encode quantum information in a way that allows for the detection and correction of errors without destroying the information itself. A prominent framework for this are topological stabilizer codes, such as the surface code, for which the logical failure rates can be systematically suppressed by increasing the size of the code if the intrinsic error rates are below some threshold value [9; 10; 11; 12; 13].

Stabilizer codes are based on a set of commutative, typically local, measurements that project an \(n\)-qubit state to a lower dimensional code space representing one or more logical qubits. Errors take the state out of the code space and are then indicated by a syndrome, corresponding to stabilizer violations. The syndrome needs to be interpreted in order to gauge whether a logical bit or phase flip may have been incurred on the logical qubit. Interpreting the syndrome, to predict the most likely logical error, requires both a decoder algorithm and, traditionally, a model of the qubit error channels. The fact that measurements may themselves be noisy, makes this interpretation additionally challenging [10; 13].

Efforts are under way to realize stabilizer codes experimentally using various qubit architectures [14; 15; 16; 17; 18; 19; 20; 21; 22; 23; 24; 25; 26; 27; 28]. In [28], code distance 3 and 5 surface codes were implemented, using 17 and 49 superconducting qubits, respectively. After initialization of the qubits, repeated stabilizer measurements are performed over a given number of cycles capped by a final round of single qubit measurements. The results are then compared with the initial state to determine whether errors have caused a logical bit- (or phase-) error. The decoder analyses the collected sets of syndrome measurements in post-processing, where the fraction of correct predictions gives a measure of the logical accuracy. The better the decoder, the higher the coherence time of the logical qubit, and in [28] a computationally costly tensor network based decoder was used to maximize the logical fidelity of the distance 5 code compared to the distance 3 code. However, with the objective of moving from running and benchmarking a quantum memory to using it for universal quantum computation, it will be necessary to do error correction both with high accuracy and in real time.

In the present work, we explore the viability of using a purely data-driven approach to decoding, based on the potential of generating large amounts of experimental data. We use a graph neural network (GNN) which is well suited for addressing this type of data. Namely, a single data point, as in [28], consists of a set of "detectors", i.e., changes in stabilizer measurements from one cycle to the next, together with a label indicating the measured logical bit- or phase-flip error. This can be represented as a labeled graph with nodes that are annotated by the information on the type of stabilizer and the space-time position of the detector, as shown in Fig
ure 1. The maximum degree of the graph can be capped based on removing edges between distant detectors, keeping only a fixed maximum number of neighboring nodes. The latter ensures that each network layer in the GNN (see Figure 2) performs a number of matrix multiplications that scales linearly with the number of nodes, i.e., linearly with the number of stabilizer measurements and the overall error rate. We have trained this decoder on simulated experimental data for the surface code using Stim [31] as well as real experimental data on the repetition code [28]. For both of these, the decoder is on par with, or outperforms, the state-of-the-art matching decoder [32], suggesting that with sufficient data and a suitable neural network architecture, model-free machine learning based decoders trained on experimental data can be competitive for future implementations of quantum error-correcting stabilizer codes.

## II Stabilizer Codes and Decoding

A stabilizer code is defined through a set of commuting operators constructed from products of Pauli operators acting on a Hilbert space of \(n\) data qubits [3]. With \(n_{S}\) independent stabilizers the Hilbert space is split into sectors of dimension \(2^{n-n_{S}}\), specified by the parity under each stabilizer. For concreteness we will consider the case \(n_{S}=n-1\), such that each of the sectors represent a single qubit degree of freedom. Each syndrome measurement is performed with the help of an ancilla qubit following a small entangling circuit with the relevant data qubits. The measured state of the ancilla qubits provide a syndrome \(S=\{s_{i},i=1,...,n_{S}\,|\in 0,1\}\), and projects the density matrix of the \(n\) qubit state into a single 2-dimensional block, a Pauli frame [33; 34]. Given uncertainties in the measurements, a number of rounds are typically performed before the information is interpreted by means of a decoder.

Defining a pair of anticommuting operators \(Z_{L}\) and \(X_{L}\) that commute with the stabilizer group, provides the logical computational space through \(Z_{L}|0\rangle_{L}=|0\rangle_{L}\) and \(|1\rangle_{L}=X_{L}|0\rangle_{L}\). Assuming a fixed pair of logical operators for a given code defines the corresponding logical states in each Pauli frame. Thus, a number of subsequent rounds of stabilizer measurements, during which the code is affected by decoherence, transforms the density matrix from the initial state \(\rho=\sum_{i,j\in\{0,1\}}\rho_{ij}|i\rangle_{L}\langle j|_{L}\) to the final state \(\rho^{\prime}=\sum_{i,j\in\{0,1\}}\rho^{\prime}_{ij}|i\rangle^{\prime}_{L} \langle j|^{\prime}_{L}\), where \(|0/1\rangle_{L}\) (\(|0/1\rangle^{\prime}_{L}\)) are the logical qubit states in the initial (final) Pauli frame. The logical error channels are approximated by

\[\rho \rightarrow\rho^{\prime}=\epsilon_{L}(\rho_{L})\] \[= (1-P)\rho+P_{X}X_{L}\rho X_{L}+P_{Z}Z_{L}\rho Z_{L}+P_{Y}Y_{L} \rho Y_{L}\,,\]

with \(Y_{L}=-iZ_{L}X_{L}\) and \(P=\sum_{i=X,Y,Z}P_{i}\). In general there may be additional non-symmetric channels (see for

Figure 1: Memory experiment on the distance \(d=5\) surface code. Data qubit initialization is followed by \(d_{t}=2\) stabilizer measurement rounds and a final data qubit measurement round. Data qubits are on the vertices of plaquettes (circles, shown in the bottom and top planes). Ancilla qubits (not shown) at the center of plaquettes provide stabilizer measurements outcomes. The detector graph has nodes corresponding to changes in stabilizers from the previous time step. (Not all edges shown.) Nodes are annotated by the type of stabilizer and the space-time coordinate. The label, here \(\lambda_{Z}=1\), corresponding to a change of \(\langle Z_{L}\rangle\), measured along the northwest edge. Also shown, bottom layer, are some example stabilizers, and the logical \(X_{L}\) (not measured).
example [19]), but we will assume that the data (as in [28]) does not resolve such channels.

The probabilities of logical error, \(P_{i}\), will be quantified by the complete set of syndrome measurements and depend on single and multi-qubit error channels as well as measurement and reset errors. It is the task of the decoder to quantify these in order to maximize the effectiveness of the error correction. Traditionally this is done through computational algorithms that use a specific error model. The framework that most decoders are based on uses independent and identically distributed symmetric noise acting on individual qubits, possibly, for circuit-level noise, complemented by two-qubit gate errors, faulty measurements and ancilla qubit reset errors. Maximum-likelihood decoders [35; 36; 37; 38; 39; 40] aim to explicitly account for all possible error configurations that are consistent with the measured syndromes, with their respective probabilities given by the assumed error model. The full set of error configurations fall in four different cosets that map to each other by the logical operators of the code, thus directly providing an estimate of the probabilities \(P_{i}\) that is limited only by the approximations involved in the calculation and the error model. Even though such decoders may be useful for benchmarking and optimizing the theoretical performance of stabilizer codes [28], they are computationally too demanding for real time operation, even for small codes.

The more standard decoders instead are based on the minimum weight perfect matching (MWPM) algorithm [41; 42; 43; 44; 45; 46]. Such a decoder aims to find the single, most likely, configuration of single qubit errors consistent with the set of measured stabilizers. Detectors are mapped to nodes of a graph with edges that are weighted by the probability of the pair of nodes. For codes where nodes appear in pairs (such as the repetition or surface code), the most likely error corresponds to pairwise matching such that the total weight of the edges is minimized. This algorithm is fast, in practice scaling approximately linearly with the size of the graph (see Figure 8). Nevertheless, it has several short-comings that limits accuracy and applicability: 1) Approximate handling of crossing edges (such as coinciding X and Z errors) means that the effective error model is oversimplified. 2) Except at very low error rates, degeneracies of less likely error configurations are ignored. 3) For models where a single error may give rise to more than two detector events, more sophisticated algorithms are needed [47; 48; 49; 50; 51; 52; 53]. These shortcomings can be partially addressed by more sophisticated approaches such as counting multiplicity or using belief propagation [54; 55; 56; 57], but often at the cost of added computational complexity. Other examples of decoder algorithms are based on decoding from small to large scale, such as cellular-automata [58; 59; 60], renormalization group [61], or union-find [49; 62]. The latter, in particular, is very efficient, but at the cost of sub-optimal performance.

### Related work

A number of different deep learning based decoder algorithms have also been formulated, based on supervised learning, reinforcement learning, and genetic neural algorithms [63; 64; 65; 66; 67; 68; 69; 70; 71; 72; 73; 74; 75; 76; 77; 78; 79; 80; 81; 82; 83]. Focusing on the works on the surface code and based on supervised learning, these can roughly be separated according to whether they primarily consider perfect stabilizers [63; 64; 65; 72; 77; 78; 81],

Figure 2: Schematic of the GNN decoder. It takes as input an annotated detector graph, c.f. Figure 1. Several layers of graph convolutional operations (following Eqn. 3) transform each node feature vector. (The empty circle shows the message passing to this particular node from neighboring nodes on the graph.) Next, a mean-pooling operation averages all the node feature vectors into a single graph embedding, which is independent of the size of the graph. Finally, the latter is passed through two separate dense networks to give two binary class predictors, corresponding to the logical \(X\) and \(Z\) labels, respectively. (For details see Appendix A.)
or include measurement noise or circuit-level noise [66, 69, 80, 82, 83], and whether they are purely data-driven [63, 65, 66, 69, 81, 83] or involve some auxiliary, model-informed, algorithm or multi-step reduction of decoding [77, 78, 80, 64, 82].

The present work is in the category, realistic (circuit-level) noise, and purely data-driven. It is distinguished primarily in that we 1) Use graph neural networks and graph structured data, and 2) Train and test the neural network decoder on real experimental data. In addition, as in several of the earlier works [66, 68, 69], we emphasize the use of a model-free, purely data-driven, approach. By using experimental stabilizer data, the approximations of traditional model-based decoder algorithms can be avoided. The fact that the real error channels at the qubit level may be asymmetric, due to amplitude damping, have long-range correlations, or involve leakage outside the computational space, is intrinsic to the data. This is also in contrast to other data-driven approaches [84, 85, 86, 21, 84] that use stabilizer data to learn the detailed Pauli channels, optimize a decoder algorithm through the edge weights of a matching decoder, or the individual qubit and measure error rates of a tensor network based decoder, as these are all constrained by a specific error model.

### Repetition code and surface code

The decoder formalism that we present in this work can be applied to any stabilizer code, requiring only a dataset of measured (or simulated) stabilizers, together with the logical outcomes. Nevertheless, to keep to the core issues of training and performance we consider only two standard scalable stabilizer codes: the repetition code and the surface code.

The repetition code is defined on a one-dimensional grid of qubits with neighboring pair-wise \(Z_{i}\otimes Z_{i+1}\) stabilizers. In the Pauli frame with all \(+1\) stabilizers, the code works are \(|0\rangle_{L}=|0\rangle^{\otimes n}\) and \(|1\rangle_{L}=|1\rangle^{\otimes n}\). Consider a logical qubit state \(|\psi\rangle=\alpha|0\rangle_{L}+\beta|1\rangle_{L}\), with complex amplitudes \(|\alpha|^{2}+|\beta|^{2}=1\). The logical bit-flip operator is given by \(X_{L}=\bigotimes_{i}X_{i}\), which sets the code distance \(d_{X}=n\). Assuming perfect stabilizer measurements and independent and identically distributed single qubit bit-flip error probabilities, decoding the repetition code is trivial. For any set of stabilizer violations, i.e., odd parity outcomes, there are only two consistent configurations of errors that map to each other by acting with \(X_{L}\). A decoder (maximum-likelihood in the case of this simple error model) would suggest the one with fewer errors. The repetition code, set up to detect bit-flip errors, is insensitive to phase flip errors, as is clear from the fact that a phase-flip (\(Z\)) error on a single qubit also gives a phase-flip error (\(\beta\rightarrow-\beta\)) on the logical qubit, corresponding to a code distance \(d_{Z}=1\). To detect and correct both bit- and phase-flip errors we need a more potent code, the most promising of which may be the surface code.

We consider the qubit-efficient "rotated" surface code [87, 88, 89] (see Figure 1), constructed from weight-4, \(Z^{\otimes 4}\) and \(X^{\otimes 4}\), stabilizers (formally stabilizer generators), with complementary weight-2 stabilizers on the boundary. On a square grid of \(d\times d\) data qubits, the \(d^{2}-1\) stabilizers give one logical qubit. We define the logical operator \(X_{L}\) as a string of \(X\)'s on the southwest edge, and a string of \(Z\)'s on the northwest edge, as shown in Figure 1. These are the two (unique up to products of stabilizers) lowest weight operators that commute with the stabilizer group, without being part of said group.

Stabilizer measurements are performed by means of entangling circuits between the data qubits and an ancilla qubit. Assuming hardware with one ancilla qubit per stabilizer, and the appropriate gate schedule, these can all be measured simultaneously, corresponding to one round of stabilizer measurements.

### Memory experiments on the surface code

To train and test our decoder we consider a real or simulated experiment, illustrated schematically in Figure 1, to benchmark a surface code as a quantum memory. The following procedure can be used for any stabilizer code:

* Initialize the individual qubits: Data qubits in a fixed or random configuration in the computational basis \(|0\rangle\) and \(|1\rangle\). Ancilla qubits in \(|0\rangle\). The initial data qubit configuration is viewed as a \(0\)'th round of measurements that initialize the \(Z\)-stabilizers. This also corresponds to an effective measurement \(\langle Z_{L}\rangle_{t=0}=\prod_{i\in Z_{L}}Z_{i}=\pm 1\). (Northwest row of qubits in Figure 1.)
* A first round, \(t=1\), of actual stabilizer measurements is performed. The Z-stabilizers are determined, up to bit-flip errors, by the \(0\)'th round. Hence, the difference between the two provides the first round of Z-detectors. The X-stabilizers have randomized outcome, projecting to an even or odd parity state over the four (or two) qubits in the Hadamard (\(|+\rangle\), \(|-\rangle\)) basis. The value of these stabilizers form the reference for subsequent error detecting measurements of the X-stabilizers. Ancilla qubits are reset to \(0\) after this and subsequent rounds.
* Subsequent rounds \(t=2,...,d_{t}\) of Z and X stabilizer measurements provide the input for corresponding detectors based on changes from the previous round.
* Finally, data qubits are measured individually in the Z-basis, which provides a final measurement, \(\langle Z_{L}\rangle_{t=d_{t}+1}\). These also provide a final round of Z-stabilizers, which, since they are provided by the actual qubit outcomes rather than by measuring an ancilla, are perfect stabilizers by definition.
The outlined experiment provides a single data point \(D=(\{V_{Z}\},\{V_{X}\},\lambda_{Z})\) consisting of set of Z-detectors \(\{V_{Z}\}\), over \(d_{t}+1\) cycles and a set of X-detectors \(\{V_{X}\}\) over \(d_{t}-1\) cycles. In addition to the stabilizer type, each detector is tagged with its space-time coordinate, \((x,y,t)\), with \(0\leq x,y\leq d\) and \(1\leq t\leq d_{t}\pm 1\) for \(Z\) and \(X\) detectors respectively. The logical label is given by

\[\lambda_{Z}=\frac{1}{2}|\langle Z_{L}\rangle_{t=0}-\langle Z_{L}\rangle_{t=d_ {t}+1}|\in\{0,1\}\,. \tag{2}\]

The probability of \(\lambda_{Z}=1\) is, according to Eqn. 1, given by \(P_{X}+P_{Y}\), and the probability of \(\lambda_{Z}=0\) by \(P_{I}+P_{Z}\), corresponding to a logical bit-flip or not.

What has been described is a "memory-Z" experiment [31], i.e., one in which we detect logical bit-flips. Qubits are initialized in the computational basis \(|0\rangle\) and \(|1\rangle\). A "memory-X" experiment prepares the qubits in the Hadamard basis, with the role of X- and Z-stabilizers reversed. Physically, in the lab, one cannot do both experiments in the same run, as \(Z_{L}\) and \(X_{L}\) do not commute. This also implies that each data point only has one of the two binary labels, \(\lambda_{Z}\) or \(\lambda_{X}\), even though there is information in the detectors about both labels. The neural network will be constructed to predict both labels for a given set of detectors, which implies that the learning framework is effectively that of semi-supervised learning, with partially labeled data. Thus, in contrast to a matching based decoder, which breaks the surface code detectors into two independent sets with a corresponding graph for each, the GNN decoder can make use of the complete information. This, in addition to the fact that it is not constrained by the limitations of the matching algorithm itself, provides a possible advantage in terms of prediction accuracy.

Additionally, some fraction of the data is incorrectly labeled. This follows from the fact that measured labels will not always be the most likely. In fact, the fraction of incorrectly labeled data corresponds to the logical failure rate that an optimal decoder would provide. For the data that we use, this fraction ranges from marginal (see Figure 4) to quite substantial (see Figure 6), depending in particular on the number of cycles that are decoded, as the logical failure rate grows exponentially with the number of cycles.

We have also assumed that there is no post-processing to remove leakage. Assuming there is some mechanism of relaxation back to the computational qubit subspace, including the last round of measurements, leakage events will be be handled automatically by the neural network decoder, based on the signature they leave in the detector data.

## III Graph neural network decoder

A graph neural network (GNN) can be viewed as a trainable message passing algorithm, where information is passed between the nodes through the edges of the graph and processed through a neural network [90, 91, 92]. The input is data in the form of a graph \(G=(V,E)\), with a set of nodes \(V=\{i\,|\,i=1,..,N\}\) and edges \(E=\{(i,j)\,|\,i\neq j\in V\}\), which is annotated by \(n\)-dimensional node feature vectors \(\vec{X}_{i}\) and edge weights (or vectors) \(e_{ij}\). The basic building blocks are the message passing graph convolutional layers, which take a graph as input and outputs an isomorphic graph with transformed feature vectors. Specifically, in this work we have used a standard graph convolution [93] where for each node the \(d_{in}\)-dimensional feature vector \(\vec{X}_{i}\) is transformed to new feature vector \(\vec{X}^{\prime}_{i}\) with dimension \(d_{out}\) according to

\[\vec{X}^{\prime}_{i}=\sigma\left(W_{1}\vec{X}_{i}+\sum_{j}e_{ij}W_{2}\vec{X}_{ i}\right)\,, \tag{3}\]

where non-existent edges are indicated by \(e_{ij}=0\). Here \(W_{1}\) and \(W_{2}\) are \(d_{out}\times d_{in}\) dimensional trainable weight matrices and \(\sigma\) is an element-wise non-linear activation function which includes a \(d_{out}\)-dimensional trainable bias vector.

For the task at hand, which is graph classification, graph convolutions are followed by a pooling layer that contracts the information to a single vector, a graph embedding, which is independent of the dimension of the graph. We use a simple mean-pooling layer \(\vec{X}^{\prime}=N^{-1}\sum_{i}\vec{X}_{i}\). For the classification we use two structurally identical, but independent, standard multi-layer feedforward networks that each end with a single node with sigmoid activation that acts as a binary classifier. The weights and biases of the complete network are trained using stochastic gradient descent with a loss function which is a sum of the binary cross entropy loss of the network output with respect to the binary labels. Since the experimental data, or simulated experimental data, only has one of the two binary labels (\(\lambda_{Z}\), \(\lambda_{X}\)) for each complete detector graph, gradients are only calculated for the provided label.

To avoid overfitting to the training data we employ two different approaches depending on the amount of available data. In using experimental data from [28], we use a two-way split into a training set and a test set. To avoid diminishing the training data further, we do not use a validation set, and instead train for a fixed number of epochs. We observe (see Figure 6) that the test accuracy does not change significantly over a large number of epochs, even though the network continues to overfits.

For the case with simulated experimental data (Figure 4), we avoid overfitting by replacing a substantial fraction (25%) of the data with new data, generated within the training cycle, after each epoch of training. This effectively emulates a much larger dataset, while keeping with the memory limits set by the available hardware. Since the training set is effectively unbounded, the number of unique detector graphs scales as \(2^{d_{i}d^{2}}\) and the network cannot overfit. Also here, a fixed test set is used to gauge the performance.
The GNN training and testing is implemented in PyTorch Geometric [94], simulated data is generated using Stim [31], and the MWPM decoding results use Pymatching [32]. The Adam optimizer is used for stochastic gradient descent, using manual learning rate decrements when the training accuracy has leveled out. Details on the training procedure can be found in Appendix A. Several other graph layers were experimented with, including graph attention for both convolutions [95] and pooling [96, 97], as well as top\({}_{k}\) pooling [98, 99]. These were found not to improve results. The width and depth of the final network was arrived at after several rounds of iterations. Naturally, we expect that larger code distances, i.e., larger graphs, will require scaling up the network size. (See also Sec. IV.4)

### Data structure

As discussed previously the data is in a form \(D=(\{V_{Z}\},\{V_{X}\},\lambda_{Z/X})\), consisting of a set of detectors \(V_{Z/X}\), specified by a space-time coordinate, together with a binary label. Based on this we construct a single graph. Each node corresponds to a detector event, and is annotated by a 5-vector (for the surface code with circuit-level noise) \(\vec{X}=(b_{1},b_{2},x,y,t)\) containing the space-time coordinate \((x,y,t)\) and two exclusive binary (one-hot encoded) labels with \(\vec{b}=(1,0)\) for an X-stabilizer and \(\vec{b}=(0,1)\) for a Z-stabilizer. (The encoding of the type of stabilizer may be superfluous, as it can be deduced from the coordinate.) We initially consider a complete graph, with edge weights given by the inverse square supremum norm between the vertices, \(e_{ij}=(\max\{|x_{i}-x_{j}|,|y_{i}-y_{j}|,|t_{i}-t_{j}|\})^{-2}\). This form of the edge weights is motivated by a naive picture of the minimum number of single data qubit measurement errors that would cause a pair of detectors. The main purpose of the weights is to give some measure of locality, in order to prune the graph. Smaller weight edges are removed, leaving only a fixed maximal node degree, which for the results presented in this work was capped at six.

## IV Results

The GNN based decoder has been implemented, trained, and tested on the surface code and the repetition code. The main focus is on using simulated or real experimental data, presented in IV.1 and IV.2, respectively. We also present some results on the surface code with perfect stabilizers, IV.3, where we are able to train the network for larger code distances.

### Surface code with circuit-level noise

We use Stim to generate data with circuit-level noise. Simulated circuits use standard settings for the surface code, containing Hadamard single qubit gates, controlled-Z (CZ) entangling gates, and measure and reset operations. All of these operations, and the idling, contain Pauli noise, scaled by an overall error rate \(p\). (See Appendix B.) Datasets of several million graphs are generated, with partial replacement after each epoch of training to avoid overfitting. Figure 3 shows test results evaluated at \(p=1.0\cdot 10^{-3}\) for decoders trained with data using an even mix of error rates \(p=\{1.0,2.0,3.0,4.0,5.0\}\cdot 10^{-3}\) and memory-Z experiments. The logical failure rate is thus approximately 50% of the true failure rate (up to correlations between failures in \(X_{L}\) and \(Z_{L}\)), but consistent with the type of data that would be experimentally accessible. (We have also tried training and testing with a mix of memory-Z and memory-X experiments, which works as well but takes longer to train to the same accuracy.) The MWPM decoder uses the information provided by the simulated error model to optimize edge weights on the decoding graph, whereas the GNN decoder uses only the data provided by the simulated measurements. Despite this, we find that with sufficient training the GNN decoder outperforms the matching decoder.

A different network is trained for each code distance \(d\) and for each number of rounds of stabilizer measurements \(d_{t}\). Figure 4 shows a representative plot of the training and validation accuracy, evaluated on the mixed error rate dataset. With an active (in-memory) dataset containing \(5\cdot 10^{6}\) and given that 25% is replaced in each epoch, 1000 epochs corresponds to a total of \(1.25\,10^{9}\) data points.

### Repetition code using experimental data

Having trained GNN based decoders on simulated experimental data in the previous section, we now turn to real experimental data. We use the public data provided together with [28]. This contains data on both the \(d=3\) and \(d=5\) surface codes as well as the \(d=25\) bit-error correcting repetition code. All datasets are of the form described in II.3, thus readily transferred to the annotated and labeled graphs used to train the GNN, as described in III.1. The datasets contain approximately \(10^{6}\) data points for the different codes, code distances, and varying number of stabilizer rounds.

Our attempts to train a GNN on the data provided for the various implementations of surface code were generally unsuccessful. While it gave good results on the training data, the logical failure rate on the test set was poor. Given the fact that on the order of \(10^{9}\) data points were used for the simulated circuit-level noise on the surface code (IV.1), it is not surprising that the significantly smaller dataset turned out to be insufficient. The net
work cannot achieve high accuracy without overfitting to the training data given the relatively small dataset.

For the repetition code, the data which is provided is of a single type, for a \(d=25\) code measured over \(d_{t}=50\) rounds. Each round thus contains the measurement of \(24\) ancilla qubits for the \(ZZ\) stabilizers of the two neighboring data qubits along a one-dimensional path. As done in [28] this data can be split up into data for smaller repetition codes, by simply restricting to stabilizers over a subset of \(d\) subsequent data qubits. In this way the dataset can be increased by a factor \(25-(d-1)\), and used to train a single GNN for each code distance. It should be noted that this is suboptimal, compared to generating the same amount of data on single distance \(d\) device, as variations in the performance of the constituent qubits and gates will be averaged out in the dataset. Nevertheless, using this scheme we successfully trained GNN decoders for short distance repetition codes, with test accuracies shown in Figure 5. Results for (what we refer to as) "Device-optimized MWPM" is taken from [28]. The GNN decoder performs almost on par with this sophisticated matching decoder for \(d=3\). As expected, the relative performance deteriorates with increased code distance. We expect that we would need more training data for larger code distance, but instead we have access to less.

As the comparison with the matching decoder that uses a device specific error model may be biased compared to using training data from different devices, as mentioned above, we also give results for an "uniformed" matching decoder with edge weights based on the 1-norm distance between space-time coordinates. It may also be noted that using MWPM corresponds to a near optimal decoder for the repetition code, at least for the case of phenomenological measurement noise where it is equivalent to bit-flip error surface code. This is in contrast to the surface code, for which MWPM is suboptimal, even in the case of perfect stabilizers. Thus, outperforming MWPM for the repetition code may be more challenging than for the surface code.

### Surface code with perfect stabilizers

To complement the results on simulated and real experimental data we have also trained the GNN decoder on the surface code with perfect stabilizers under depolarizing noise. The same network (see Appendix A) is used as for circuit-level noise, but trained at \(p=[0.01,0.05,0.1,0.15]\). Results up to code distance \(d=15\) are shown in Figure 7 and found to significantly outperform MWPM. We also compare to a tensor network based [100] maximum likelihood decoder (MLD), showing that for code distance \(d=5\) the GNN decoder has converged to the level of being an approximate MLD.

We do not attempt to derive any threshold for the GNN decoder. Given a sufficiently expressive network we expect that the decoder would eventually converge to a maximum likelihood decoder, but in practice the accuracy is limited by the training time. It gets progressively more difficult to converge the training for larger code distances, which means that any th

Figure 4: GNN training and test accuracy versus number of training epochs for circuit-level noise, comparing training using a fixed dataset to replacing \(25\%\) of the data after each epoch. We see that the latter setup is crucial to avoid overfitting, allowing the test accuracy to closely follow the training accuracy. Code distance \(d=5\), with \(d_{t}=5\) cycles, dataset size \(5\cdot 10^{6}\) samples using an error rate randomly selected from \(p=[0.001,0.002,...,0.005]\). The test set is a fixed dataset of the same type containing \(5\cdot 10^{4}\) data points. (Sharp kink at \(250\) epochs correspond to a decrement of the learning rate.)

Figure 3: Decoding simulated experimental data [31], with circuit-level noise, on the rotated surface code with code distance \(d\). Logical failure rate versus number of rounds of stabilizer measurements \(d_{t}\). Comparing GNN decoder trained on the detector data with MWPM decoder [32] that has full information of the data-generating error model. Each data point is evaluated over \(10^{8}\) samples (\(10^{7}\) for \(d<7\)), memory-Z, at a reference error rate of \(p=1\cdot 10^{-3}\).
a function of the training time versus code distance. In fact, in principle, since the threshold is a \(d\to\infty\) quantity, we would not expect that a supervised learning algorithm can give a proper threshold if trained separately for each code distance. Nevertheless, through GNN's it is quite natural to use the same network to decode any distance code, as the data objects (detector graphs) have the same structure. We have investigated training the same network for different code distances and different number of rounds. This shows some promise, but so far does not achieve accuracy levels that can match MWPM.

### Scalability

We are limited to relatively small codes in this work. For the repetition code using experimental data, it is quite clear that main limitation to scaling up the code distance is the size of the available dataset. For the surface code using simulated data it is challenging to increase the code distance while still surpassing MWPM. As the logical failure rates decrease exponentially with code distance, the test accuracy of the supervised training needs to follow. One way to counter this is to increase the number of stabilizer cycles, \(d_{t}\), but this also increases the graph size, making the training more challenging from the perspective of increased memory requirements as well as the increased complexity of the data.

Nevertheless, it is interesting to explore the intrinsic scalability of the algorithm, by quantifying how the decoding time using a fixed size GNN scales with the code size. Here we present results on the decoding time per syndrome for the surface code, as a function of code volume \(d^{2}d_{t}\), at fixed error rate. The network architecture is the same as used for all the results in this paper, as described in Appendix A. In line with expectations, the GNN inference scales approximately linearly with the code volume, i.e. average graph size, \(T\sim d^{2}d_{t}\) (with

Figure 5: Decoding experimental data [28] on the repetition code with code distance \(d\), over 50 rounds of stabilizer measurements. Comparing GNN decoder, using a dataset containing \((26-d)\cdot 5\cdot 10^{7}\) graphs, with a MWPM decoder with “device-optimized” edge weights ([28]) and a simple model-free MWPM decoder with 1-norm edge weights. The training-test split of the dataset is 99 to 1, and the logical failure rate is mapped to an error rate per round. Results for two different random training-test splits are shown.

Figure 6: Training curves for the GNN decoder on repetition code, following Figure 5. Each epoch trains through the whole training set, which eventually leads to overfitting, where the training accuracy starts to significantly surpass the test accuracy. To maximize the amount of training data, no validation set was used. No early stopping was implemented in order to avoid optimizing results to the test set.

Figure 7: Decoding the rotated surface code with perfect stabilizers and code distance \(d\). Logical failure rate versus error rate \(p\), for depolarizing noise, evaluated over failures with respect to both \(X_{L}(\lambda_{X})\) and \(Z_{L}(\lambda_{Z})\). Comparing the GNN decoder with MWPM decoder that has full information of the data-generating error model. Each data point is evaluated over \(10^{5}\) data points. Dashed lines is the accuracy using a matrix product state (MPS) decoder [100] at code distances 5 and 7, evaluated over \(10^{4}\) data points.
\(d_{t}=1\) for perfect stabilizers). The number of matrix operations per graph convolutional layer, following Equation 3, is proportional to the number of nodes times the number of edges. The number of layers is fixed, multiplying this by a constant factor. The feature vector pooling is proportional to the number of nodes, whereas the subsequent dense network classifiers are independent of the graph size. We find that inference scales slightly better than the highly optimized matching decoder. However, several caveats are in order. 1) The size of the GNN is fixed. Larger code sizes may eventually require scaled up networks, unless the error rate is scaled down accordingly 2) The network has not been trained on code distances larger than \(d=15\) (2D). It is only a test of the decoding time, not the accuracy. 3) For GNN inference, the data is batched in order to eliminate a fixed loading time to the GPU. Treating batched data doesn't seem viable for real time decoding. Similarly our graph construction algorithm is slower, scaling quadratically with code volume, and this time has also been removed to get decoding time per graph. These are both issues that most likely can be improved significantly with more hardware efficient algorithms, and, in the longer term, special purpose hardware.

## V Conclusion and Outlook

In this paper we develop a model-free, data-driven, approach to decoding quantum error correcting stabilizer codes, using graph neural networks. A real or simulated memory experiment is represented as a single detector graph, with annotated nodes corresponding to the type of stabilizer and its space-time coordinate, and labeled by the measured logical operation on the encoded qubit. The maximal node degree is capped by cropping edges between distant nodes. The data is used to train a convolutional GNN for graph classification, with classes corresponding to logical Pauli operations, and used for decoding. We show that we can use real and simulated experimental data,for the repetition code and surface code respectively, to train a decoder with logical failure rates on par with minimum weight perfect matching, despite the latter having detailed information about the underlying real or simulated error channels. The use of a graph structure provides an efficient way to store and process the syndrome data. To train the GNN requires significant amounts of training data, but as shown in the case of simulated experiments, data can be produced in parallel with training. Network inference, i.e., using the network as a decoder, is fast, scaling approximately linearly with the space-time dimension of the code.

As an extension of this work there are several interesting possibilities to explore. One example is to use a GNN for edge weight generation within a hybrid algorithm with a matching decoder (similarly to [21]). This would depart from the pure data-driven approach pursued in this paper, with peak performance limited by the matching decoder, but with the potential advantage of requiring less data to train. An alternative to this, to potentially improve performance and lower data requirements, is to use device specific input into edge weights, or encode soft information on measurement fidelities into edge or node features.

Going beyond the most standard codes considered in this paper, we expect that any error correcting code for which labeled detector data can be generated can also be decoded with a GNN. This includes Clifford-deformed stabilizer codes [101, 102, 103, 104, 105], color codes [106, 107] or hexagonal stabilizer codes [108, 109, 110, 111], where syndrome defects are not created in pairs, but potentially also Floquet type codes [112, 113]. In addition, heterogeneous and correlated noise models [114, 115] would also be interesting to explore, where in particular the latter is difficult to handle with most standard decoders.

The software code for the project can be found at [116].

###### Acknowledgements.

 We acknowledge financial support from the Knut and Alice Wallenberg Foundation through the Wallenberg Centre for Quantum Technology (WACQT). Computations were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS) and the Swedish National Infrastructure for Computing (SNIC) at Chalmers Centre for Computational Science and Engineering (C3SE), partially funded by the Swedish Research Council through grant agreements no. 2022-06725 and no. 2018-05973. We thank

Figure 8: Scaling of average decoding time per syndrome versus code volume \(d^{2}d_{t}\) for GNN and MWPM. Perfect stabilizer measurements sampled at \(p=0.05\) (labeled 2D) and circuit level noise sampled at \(p=10^{-3}\) (labeled 3D). Dotted lines show a regression for \(d\geq 31\) and \(d\geq 11\) for 2D and 3D, respectively, according to the ansatz: \(T=C\cdot(d^{2}\cdot d_{t})^{\alpha}\), with the GNN showing sub-linear scaling. The average graph construction time per sample is shown in green.
Viktor Rehnberg and Hampus Linander for technical support.

## Appendix A GNN architecture and training

Figure 9 displays the architecture of the GNN decoder. The node features are sent through 7 subsequent graph convolutional layers (Equation 3). The node features are passed through a rectified linear unit (ReLU) activation function (which corresponds to chopping negative values) after each layer. After the graph convolutional layers, the node features from all nodes are pooled into one high-dimensional vector by computing the mean across all nodes. This vector is then cloned and sent to two identical fully connected neural networks. Both heads consist of 4 dense layers which map the pooled node feature vector down to one real-valued number which is output in the range 0 to 1 through a sigmoid function. The input and output dimension \(d_{in}\) and \(d_{out}\) of the graph convolutional and dense layers can be found in Table 1.

Networks are trained on NVIDIA Tesla A100 HGX GPU's using the python multiprocessing module to generate data in parallel on a CPU. For gradient descent, samples are batched in batches of size \(10^{3}\). The learning rate is set to \(10^{-4}\) and decreased manually to \(10^{-5}\), whenever the validation accuracy reached a plateau. An example of a training history for \(d=5\) and varying number of surface code cycles \(d_{t}\) is shown in Figure 10. For this example, with \(d_{t}=5\), 100 epochs of training takes approximately 10 hours. The code is available at [116].

## Appendix B Stabilizer circuits and error model for circuit-level noise

Quantum circuits for weight-four \(Z\)- (\(X\)-) stabilizers of the surface code are displayed in Figure 11 (12). The gate set used for the stabilizer measurements consists of the Hadamard gate (\(H\)), and the \(CNOT\) gate. Under

\begin{table}
\begin{tabular}{|c|c|c|} \hline \hline Layer & \(d_{in}\) & \(d_{out}\) \\ \hline GraphConv\({}_{1}\) & 5 & 32 \\ GraphConv\({}_{2}\) & 32 & 128 \\ GraphConv\({}_{3}\) & 128 & 256 \\ GraphConv\({}_{4}\) & 256 & 512 \\ GraphConv\({}_{5}\) & 512 & 512 \\ GraphConv\({}_{6}\) & 512 & 256 \\ GraphConv\({}_{7}\) & 256 & 256 \\ Dense\({}_{1}\) X & 256 & 128 \\ Dense\({}_{2}\) X & 128 & 181 \\ Dense\({}_{3}\) X & 128 & 32 \\ Dense\({}_{4}\) X & 32 & 1 \\ Dense\({}_{1}\) Z & 256 & 128 \\ Dense\({}_{2}\) Z & 128 & 181 \\ Dense\({}_{3}\) Z & 128 & 32 \\ Dense\({}_{4}\) Z & 32 & 1 \\ \hline \end{tabular}
\end{table}
Table 1: Overview over the input and output dimension of the graph convolutional and dense layers of the GNN decoder.

Figure 10: Example of training and validation accuracy versus the number of epochs of training for code distance 5. One epoch corresponds to training with a dataset containing \(5\cdot 10^{6}\) detector graphs of different error rates, as in Figure 4. The test set is a fixed dataset of the same type containing \(5\cdot 10^{4}\) data points. After each epoch the oldest 25% of the dataset is replaced with new data to avoid overfitting to the training data. Kink at 300 epochs correspond to decrement of the learning rate, whereas the spikes are due to network fluctuations.

Figure 9: Schematic of the GNN architecture, with details in Table 1. The same architecture is used for all the results, except that for the repetition code there is only one output head. Also the input dimension is two (2D space-time coordinate) for the repetition code and four (two types of stabilizers, and 2D spatial coordinate) for the surface code with perfect stabilizers.
circuit-level noise, single-qubit depolarizing noise gate \(D_{p}\) (which applies gate \(\sigma_{i},i\in\{X,Y,Z\}\) where any of the gates is applied with probability \(p/3\), and \(I\) with probability \(1-p\)) acts on the data qubits before each stabilizer measurement cycle and on each target qubit after single-qubit gates. Two-qubit depolarizing noise gates (which apply gate \(\sigma_{i}\sigma_{j},i,j\in\{I,X,Y,Z\}\), where \(II\) is acted on with probability \(1-p\), and the rest with probability \(p/15\)) act on the two qubits involved after every two-qubit gate. Furthermore, each qubit suffers from reset- and measurement-error with probability \(p\), displayed by operators \(X_{p}\) when measuring and resetting in the computational basis.
