The author wants to thank Ben Criger, Barbara Terhal, Thomas O'Brien for useful discussion.
The implementation of minimum-weight perfect matching algorithm, including the one used in [appendix:varyingperror]Appendix , is provided by Christophe Vuillot, which uses the backend from either Blossom V [1]} or NetworkX [2]}.
The author acknowledge support through the ERC Consolidator Grant No. 682726.


Implementation of Belief Propagation Algorithm
Belief propagation is a heuristic procedure for computing marginal probabilities of graphical models.
We choose to use a slightly different belief propagation implementation compared to [3]}, as ours seems to be more natural for the bit-flip noise model.
We divide the lattice into \(2\times 2\)  unit cells.
Let \(G\)  be a bipartite graph, where one part corresponds to unit cells, and another part to coarse-grained edges.
Two vertices in \(G\)  is connected when the coarse-grained edge is adjacent to the unit cell.
This later decides how the messages flow in the graph.
However, we assign two variables \(\lbrace x(e_i),x(e_j)\rbrace \)  to each vertex corresponding a coarse-grained edge which \(e_i\)  and \(e_j\)  form.
In this section, the symbol \(e\)  or \(e_i\)  will denote original edges of the lattice (i.e. not coarse-grained).
We define \(E\)  to be the set of all edges \(e\) , \(E_{cg} \subset E\)  to be the set of \(e\)  which are components of coarse-grained edges (i.e. red edges in fig:lattice), and \(\bar{E}_{cg} = E \setminus E_{cg}\) .
Given a syndrome \(S\) , the unnormalized probability of an error configuration \(\vec{x}\equiv \lbrace x(e)\rbrace _{e\in E}\)  can be written as
\(p(\vec{x}) = g(S, \vec{x})\prod _{e\in E} p_e(x(e)) ,\) 
where \(g(S, \vec{x}) = 1\)  if \(\vec{x}\)  has syndrome \(S\) , and otherwise \(g(S, \vec{x}) = 0\) .
It is obvious \(g(S, \vec{x})\)  can be factorized to local terms.
Thus, the marginal distribution for \(e\in E_{cg}\)  can then be factorized according to \(G\)  as
\(\sum _{\lbrace x(e), e\in \bar{E}_{cg}\rbrace }p(\vec{x}) = \prod _c f_c(\lbrace x(e)\rbrace _c) ,\) 
where the product is taken over all unit cells \(c\) , and \(\lbrace x(e)\rbrace _c\)  are the set of \(x(e)\)  such that \(e\in E_{cg}\)  is adjacent to \(c\) .
We can then apply the standard belief propagation to the graph \(G\) .
Without further explanation, we choose to use the following rule.
A unit cell \(c_k\)  sends to each of its adjacent cell \(c_n\)  a message containing 4 real numbers
\(\lbrace m_{c_k, c_n}(x(e_i), x(e_j))\rbrace  \quad \text{for} \ x(e_i),x(e_j) = 0,1 ,\) 
where \(e_i\)  and \(e_j\)  form the coarse-grained edge between \(c_k\)  and \(c_n\) .
When we already fix an error configuration \(\vec{x}\) , we may use the simplified notations
\(p_e \equiv p_e(x(e)), \qquad m_{c_k, c_n} \equiv m_{c_k, c_n}(x(e_i), x(e_j)).\) 
To compute an out-going message from cell \(c\) , we take messages from the other three directions, and consider them to be the probability of error configuration on respective edges.
We then sum over all error configurations in the cell, which give the correct syndrome of the 4 plaquette stabilizer checks.
<FIGURE>More concretely, we define \(\vec{x}_c\)  to be \(\vec{x}\)  restricted to edges in \(c\)  (i.e. all edges in fig:bpcell), and \(g^{\prime }(S, \vec{x}_c)\)  checks whether \(\vec{x}_c\)  is compatible with \(S\)  similar to \(g\) .
We assume we want to send the messages from \(c\)  to \(c_n\) , while the incoming messages are from cells \(I = \lbrace c_k\rbrace \) .
Then we have
\(m_{c, c_n} = \sum _{\vec{x}_c} g^{\prime }(S, \vec{x}_c)\prod _{c_k\in I} m_{c_k, c} \prod p_{e_i}.\) 
For the last term \(\prod p_{e_i}\) , the product is taken over the blue edges in fig:bpcell, assuming \(c_n\)  is on the right of \(c\) .
In the end of the message passing, we can compute the marginal probability by
\(P(x(e_i),x(e_j)) = m_{c_k, c_n}m_{c_n,c_k} / \left(p_{e_i}p_{e_j}\right) ,\) 
where \(e_i\)  and \(e_j\)  are the edges between \(c_n\)  and \(c_k\) .
From the joint distribution \(P(x(e_i),x(e_j))\)  we can compute the distribution \(P\left(x(e_i)+x(e_j)\mod {2}\right)\) .
It is not hard to see that the above message passing rules will lead to the correct marginal probability when the underlying graph is a tree
(note this is not the case for the square lattices we are considering).
To generate training data for neural networks, we do 7 rounds of message passing defined above.
The key differences between our implementation and the one in [4]} is

Ours utilizes all four stabilizer checks in each unit cell while in [4]} only three are used.

Each message contains 4 real numbers in our implementation while only 1 in [4]}.


Introduction to Neural Networks
A neural network, at the highest abstraction, can just be viewed as a black-box function \(f_{\text{nn}}(x,\vec{w})\)  with many parameters \(\vec{w}\)  to be tuned.
We want \(f\)  to describe the input-output relation presented in a dataset \(D=\lbrace (x_i, y_i)\rbrace \) .
To do thisIn this work, we choose to not study the effect of overfitting, as we have the ability to generate infinite data., we choose a (smooth) loss function \(\mathcal {L}\) , and then we do the minimization
\(\min _{\vec{w}}\sum _i \mathcal {L}(f(x_i, \vec{w}), y_i) .\) 
One important requirement is that \(f\)  is (almost-everywhere) differentiable with respect to \(\vec{w}\) .
This allows us to train the network with gradient descent, for which a good introduction can be found in [7]}.
In general, we can expect gradient descent will take us to a local minimum or some region with very small gradients.
This is the advantage of “end-to-end” training compared to human-written heuristic algorithms, as the latter are unlikely to be a local optimum (assuming we can add real number parameters to those heuristic algorithms).
A common loss function for classification problems is cross-entropy.
Assume we have a dataset \(D=\lbrace (x_i, y_i)\rbrace \)  where \(y_i\in \lbrace 0,1\rbrace \) , and the neural network output \(y_i^{\prime }\)  which tries to approximate \(\text{Prob}(y_i=1)\) , the cross-entropy loss function is then calculated as following:
\(-\sum _i \left( y_i \log y_i^{\prime } + \left(1-y_i\right) \log \left(1-y_i^{\prime }\right)\right).\) 
Note that when we use the notation \(\text{Prob}(y_i=1)\) , we implicitly assume \(D\)  is obtained by sampling from an underlying probability distribution.
More concretely, most neural networks consist of many layers.
In this paper, the two relevant types of layers are the dense and convolution layer.
Dense layers (a.k.a fully-connected layers) have the form \(g(A\vec{x}+\vec{b})\) , where \(g\)  is some non-linear function applied entrywise, and the matrix \(A\) , vector \(\vec{b}\)  are the trainable parameters.
Assuming \(A\)  has a shape of \(n\times m\) , we will say the output dimension of the dense layer is \(n\) .
One convolution layer, as the name suggests, contains a collection of discrete convolutions.
For this paper, the input to the layer resides on a 2-dimensional lattice of size \(l^2\)  with periodic boundary condition.
On each lattice site, there is a \(d\) -dimensional input vector \(x_{\vec{u}} \in \mathbb {R}^d\) , where the subscript \(\vec{u} \in \mathbb {Z}_l^2\)  (we use \(\mathbb {Z}_l\)  to denote integer in range \([0,l-1]\) ).
We define the kernel to be a tensor \(K_{\vec{u},i}\) , where \(\vec{u} \in \mathbb {Z}_n^2\)  and \(i \in \mathbb {Z}_d\) . With a slight abuse of notation, we will say such a kernel has size \(n^2\) . The convolution is then
\(y_{\vec{v}} = \sum _{\vec{u} \in \mathbb {Z}_n^2} \sum _{i \in \mathbb {Z}_d} x_{\vec{v}-\vec{u}, i} K_{\vec{u}, i},\) 
where \(x_{\vec{v}-\vec{u},i}\)  is the \(i\) th element of \(x_{\vec{v}-\vec{u}}\) , and \(\vec{v}-\vec{u}\)  is calculated module \(l\)  because of the periodic boundary condition.
We will have a collection of kernels \(\lbrace K_{\vec{u}, i,j}\rbrace _j\)  for one convolution layer, which means we also have a collection of outputs \(\lbrace y_{\vec{v},j}\rbrace _j\) .
The cardinality of \(\lbrace K_{\vec{u}, i,j}\rbrace _j\)  is conventionally called the number of filters.
After eq:convolution, we can also apply a non-linear function \(g\)  entrywise if needed.
Before concluding this section, let us make one clarification.
In this paper, sometimes we only train part of the network, e.g. for blue circle 1 in subfig:completenetwork we only train the dense layers.
Assuming \(\vec{w}_1\)  are the parameters in the part of the network we want to train and \(\vec{w}_2\)  are the rest, then we are doing the optimization
\(\min _{\vec{w}_1}\sum _i \mathcal {L}(f(x_i, \vec{w}_1, \vec{w}_2), y_i)\) 
by using some gradient descent optimizer.

Comparison to Simpler Approachers
In this section, we will show the performance of the neural net decoders when trained with simpler approaches (more precisely, approaches with less human involvement and prior knowledge of toric code decoding), and provide some reasoning if possible.
The neural nets will be the same as the ones we used in the main text, except that they do not contain the “removing complexity” step in the post-processing.
We will see in general the performance gets much worse, especially when the lattice size grows large.
However, this does not mean these simpler approaches will always fail.
It just implies that a large amount of training time / human involvement is needed, which could make them impossible in practice.
The simplest approach is to train the whole network with input-output pairs (syndrome, logical correction).
During limited attempts, this approach does not produce decoders much better than random guess for large toric codes.
A hand-waving explanation is the following.
It is fair to assume a lot of parity functions need to be evaluated during the decoding process.
It is known that the parity is not an easy function for neural nets to compute [8]}, and one good way to approximate it is to increase the depth of the network.
So let us assume each renormalization stage needs 5 layers.
This means to decode \(L=32\)  toric code, the network will have 25 layers, which exceeds the range where neural nets can be reliably trained.
The problem of too many layers can be alleviated if we can pre-train the earlier layers of the network.
A similar strategy was used in training neural nets for computer vision problems [9]}.
For the bit-flip noise model, we can pre-train the earlier layers to mimic the renormalization group decoder.
Recall in sec:toriccodeandrndecoder, we mentioned that for the renormalization group decoder, the output corresponding to a coarse-grained edge \(e\)  is \(p_e\) .
Since we generate syndromes by first sampling \(x(e)\)  for all (not coarse-grained) edges \(e\) , we also have the ability to generate pairs (syndrome, \(\lbrace x(e)\rbrace \)  for coarse-grained edges \(e\) ) for training.
Although \(\lbrace x(e)\rbrace \)  in the above pairs are binary numbers, with the cross-entropy as the cost function, in theory the output will converge to \(p_e(1)\) .
The pre-training is done one renormalization block at a time.
More concretely, with the network we are using in the main text, we will train the output of the 12th layer with the training target of first renormalization block, and the output of the 24th layer with the second block, etc.
We can try this method on \(L=32\)  toric code and bit-flip error rate \(0.08\) .
It does not work well, as we can see in tab:crossentropyafterblocks.
<TABLE>These numbers are not very accurate and coming from a single training instance.
However, the author has observed the same trend many times that the loss and accuracy slowly degrade in the process of renormalization, even though the error rate is way below the theoretical threshold.
To diagnose the reason, we first notice that the first RN block actually performs reasonably well.
This suggests that the optimizer is capable of training each RN block alone.
Assume this is indeed the case, the degrading of performance is likely caused by the following two reasons.
First, later in the renormalization process, if we look at the coarse-grained syndrome or \(p_e\)  alone, they behave more and more like white noise.
While it is possible for the neural nets to do the same post-processing described in subsec:designnetwork, a few layers of the network will be occupied by this.
Therefore, the natural solution is to implement the post-processing ourselves.
By doing this, we suspect it is possible to reach a threshold of \(8\%\) , but apparently \(8\%\)  is still not good enough.
The second reason is related to the convergence of \(p_e\) .
When below the threshold, we will encounter very often that \(p_e\)  is very close to 0 or 1.
For example, if \(x(e)=1\)  vs \(x(e)=0\)  corresponding to a weight 4 vs weight 1 local configuration, then we will have \(p_e \approx p_0^3\) , where \(p_0\)  is the initial error rate.
This will become more prominent later in the renormalization process, as \(p_e\)  from the previous renormalization stage become the error rate in the next stage.
It is important to know how close to 0 (or 1) \(p_e\)  is on the logarithmic scale.
Otherwise, in the later stage, the information will not be accurate enough to deduce configuration close to the minimum weight of errors.
This poses the following requirements:

When we pass \(p_e\)  to some intermediate layer of a neural net, it should be able to distinguish between small \(p_e\) .
However, recall that each layer does the computation \(f(Ax + b)\) , where \(f\)  has a bounded derivative.
Thus, to distinguish a set of small \(\lbrace p_e\rbrace \) , we need \(\Vert  A \Vert  \sim 1/\min \lbrace p_e\rbrace \) .
This will either not be achieved by training or cause instability of the network.
Another issue related to the minuscule nature of \(p_e\)  is the cross-entropy loss function does not provide enough motivation for \(p_e\)  to converge to the target value \(q\)  in log-scale.
More accurately, the derivative of the cross-entropy scales like \(O(|p_e - q|)\)  when \(p_e \approx q\) , which will be too small before the convergence in log-scale.
A natural solution is we replace the appearance of \(p_e\)  with \(\log p_e -\log (1-p_e)\) .

Even with a good representation of \(p_e\) , we shall still be very cautious about the training, as we are trying to estimate very small \(p_e\)  from sampling.
In the end, we decide to implement a belief propagation routine and use the input-output pair from the routine to train the network.
The advantage is that belief propagation directly outputs the probability, which should be reasonably accurate in the logarithmic scale.
Therefore, we can get a much stable training process.


Technical Details
The objective of this section is to describe some technical details for people who do not plan to read the source code.
For the majority of network, we use leaky ReLUs [10]} as the activation function, which has the form
\(y=x \quad \text{if} \ x>0; \quad y=0.2x \quad \text{if} \ x\le 0\) 
Apart from the last layer of each renormalization block and the last dense layer, the number of filters in each convolution layer is 200, and the output dimension of each dense layer is 50.
For training a belief propagation network, we generate a dataset of size 80000.
The dataset consists of the input and output of the belief propagation algorithm described in [appendix:bp]Appendix  when applied to \(d=16\)  toric code.
The optimizer we use is the ADAM [11]}.
The learning rate parameter of the optimizer is set to \(7\times 10^{-4}\)  for training belief propagation network, \(10^{-3}\)  for training dense layers, \(7\times 10^{-5}\)  for global training of \(L=16,32\)  lattice, and \(7\times 10^{-6}\)  for \(L=64\) .
The batch size for training is 50.
The training of the dense layers uses around 1000 batches.
For the strong policy, the global training uses 3000 batches regardless of the code distance.
To see the potential of our neural decoder at \(d=64\) , we also did a longer training and compared it to MWPM in fig:comparisontoMWPM.
In total, it is trained using 18000 batches.
The first 12000 batches are trained on physical error rate \(p=0.09\) , and the last 6000 batches are on \(p=0.095\) .
The reason of switching to a higher error rate for late-stage training is that the accuracy at \(p=0.09\)  is too close to 1 for effective training.

Spatially Varying Error Rates
A natural use of the neural decoder is to train it with experimental data.
Naturally, the noise models in the experiments will not be translational invariant.
There are two simple ways to reconcile this fact with the translational-invariance of convolutional neural nets:

Allow the first few layers of the network to be non-translational invariant.

Introduce site-dependent trainable variables to the networks.

In this section, we will consider the error model that has varying bit-flip error rates across the lattice.
For this, we can use the second approach, where the site-dependent variables can in principle represents the varying error rates.
In fact, recall that the neural decoder has an input shape \((L,L,3)\) , which contains \(2L^2\)  numbers that are originally error rates feeding into belief propagation.
Thus we can simply feed the site-dependent variables into the neural decoder and then train them.
However, there is still a complication regarding the starting point of this training.
A natural choice is to start with the trained neural decoder for uniform error rate and only train the site-dependent variables.
With this route, there is a risk that if previously we trained the neural decoder under uniform error rate for long enough, the neural decoder could learn to ignore the error rate inputs as they are constant.
In this case, only training the site-dependent variables can fail.
Another route we can take is that we also reinitialize the first renormalization block.
More accurately, we do the following:

We start with the trained neural decoder for uniform error rate.

We reverse the belief propagation network in the first renormalization block to the pre-trained weights. In other words, it now again approximates belief propagation.

We then train the site-dependent variables and the first renormalization block together.

We test this scheme with a distance 16 toric code.
For each qubit, there is a \(50\%\)  chance it has bit-flip rate \(0.16\) , and a \(50\%\)  chance of rate 0.
In other words, only around half of the qubits are noisy.
With this new error model, we train the site-dependent variables and the first renormalization block for 4500 batches, where each batch contains 50 training data.
The learning rate of the ADAM optimizer is set to \(2\times 10^{-4}\) .
By performing this training, we increase the logical accuracy from \(0.967\)  to \(0.993\) , where \(0.967\)  corresponds to the decoder trained on uniform bit-flip error model.
The accuracy is each evaluated by using \(10^5\)  (syndrome, logical correction) pair.
To provide some comparison, we run the minimum-weight perfect matching (MWPM) algorithm for the same error model.
Without providing error rate information, the MWPM algorithm assigns an equal weight to each qubit, and have a logical accuracy of \(0.975\) .
If we provide the perfect information about the error rates, we can assign weight 1 to noisy qubits, and weight 100 to the noiseless ones.
For each pair of violated parity checks, we can select the path with the minimum total weight between them, and use this weight for the MWPM algorithm.
With this choice, the logical accuracy rises to 1, e.g. \(100\%\)  success.
These accuracies are each evaluated by \(10^4\)  (syndrome, logical correction) pair.
Based on the thoughts above, it is likely better to not start with the neural decoder trained on uniform error rate.
Instead, we can train a neural decoder with training data that has varying error rates, but otherwise the same procedure as depicted in subfig:completenetwork.
This way, the first renormalization block will not learn to ignore the error rate inputs.
