Sedain & Kacharov IAU Symposium 379

17 2023 10.1017/xxxxx

Proceedings of IAU Symposium 379

# Testing Jeans dynamical models with prolate rotation on a cosmologically
simulated dwarf galaxy

Amrit Sedain1,2 Nikolay Kacharov1 1 Leibniz Institute for Astrophysics, An der
Sternwarte 16, 14482 Potsdam, Germany
2 Institute of Physics and Astronomy, Karl-Liebknecht Str.24/25 14476 Potsdam,
Germany

###### Abstract

Prolate rotation is characterized by a significant stellar rotation around a
galaxy’s major axis, which contrasts with the more common oblate rotation.
Prolate rotation is thought to be due to major mergers and thus studies of
prolate-rotating systems can help us better understand the hierarchical
process of galaxy evolution. Dynamical studies of such galaxies are important
to find their gravitational potential profile, total mass, and dark matter
fraction. Recently, it has been shown in a cosmological simulation that it is
possible to form a prolate-rotating dwarf galaxy following a dwarf-dwarf
merger event. The simulation also shows that the unusual prolate rotation can
be time enduring. In this particular example, the galaxy continued to rotate
around its major axis for at least $7.4$ Gyr (from the merger event until the
end of the simulation). In this project, we use mock observations of the
hydro-dynamically simulated prolate-rotating dwarf galaxy to fit various
stages of its evolution with Jeans dynamical models. The Jeans models
successfully fit the early oblate state before the major merger event, and
also the late prolate stages of the simulated galaxy, recovering its mass
distribution, velocity dispersion, and rotation profile. We also ran a
prolate-rotating N-body simulation with similar properties to the
cosmologically simulated galaxy, which gradually loses its angular momentum on
a short time scale $\sim 100$ Myr. More tests are needed to understand why
prolate rotation is time enduring in the cosmological simulation, but not in a
simple N-body simulation.

###### keywords:

Galaxy dynamics, Prolate rotation, Cosmological simulations

††editors: P. Bonifacio, M.-R. Cioni & F. Hammer, eds.

##  1 Introduction

Dwarf galaxies are important building blocks in galaxy hierarchical formation
and evolution theories, so it is essential to understand their properties.
There are distinct types of dwarf galaxies, depending on mass and gas content
- mainly gas-rich dwarf irregulars and gas-poor dwarf spheroidals, as well as
transition types. The majority of them are dark matter (DM) dominated.
However, their low stellar densities and shallow potential wells make them
very sensitive to perturbations from stellar feedback, external interactions,
etc.

Dynamical studies of dwarf galaxies can reveal their total mass and the mass
of their sub-components (stars, gas, DM content), as well as their internal
structure and kinematical properties, which is not possible from photometry
alone.

On the other hand, cosmological simulations allow us to study the formation
and evolution of galaxies at all scales and test real-time observations. By
comparing observations and the simulated results we can better understand the
physical processes in our Universe. Zoom-in simulations are used to study
smaller-scale processes in individual galaxies and galaxy clusters, such as
interactions, gas cooling, star formation, and feedback from supernovae and
black holes.

While most galaxies have a certain amount of angular momentum, which flattens
them in the direction of the rotation axis, there are rare exceptions of
galaxies that rotate around their major axis, known as prolate rotation. In
the Local Group, we know of only two such cases, the And II dwarf (Amorisco et
al., 2014) and the Phoenix dwarf (Kacharov et al., 2017). In both studies,
they strongly suggest that a major merger at some point in the galaxy’s
evolution causes a change in the rotation axis. Similarly, Cardona-Barrero et
al. (2021) discovered the presence of prolate rotation in a cosmologically
simulated dwarf galaxy, using the gear code, based on gadget-2 (Springel,
2005; Revaz & Jablonka, 2012). The culprit, causing the flip of the rotation
axis is also a major dwarf-dwarf merger with a DM mass ratio of about $1:5$.

In this study, we aim to test how well Jeans dynamical models can reproduce
the mass distribution and kinematic properties of the simulated prolate-
rotating dwarf. If successful, such dynamical models can then be applied to
real observational data.

##  2 Jeans Anisotropic Multi-Gaussian Expansion (JAM) on the cosmological
data

Jeans equations are stellar hydrodynamic equations, derived from the
collisionless Boltzmann equation, that are used to study the dynamics of
equilibrium systems. One of the advantages of the Jeans equations is that
instead of dealing with the phase-space distribution function, we use velocity
moments of the line-of-sight velocity distribution (LOSVD), which is directly
linked to the gravitational potential and density of the system. This method
has been extensively used to study spherical and axisymmetric stellar systems
under the assumption of quasi-equilibrium (Cappellari, 2016).

In this work, we utilise the JAM code by Cappellari (2020) to solve the Jeans
equations for the cosmologically simulated prolate rotating dwarf galaxy,
discovered by Cardona-Barrero et al. (2021). We obtain best-fit Jeans
dynamical models at different stages during its evolution in the cosmological
simulation, which predicts the dwarf’s mass distribution and compares the
model results with the known properties of the simulated dwarf. The primary
objective is to investigate whether the Jeans models can accurately describe
the dynamical properties of the prolate system in a quasi-equilibrium
approximation since in the cosmological simulation the prolate rotation
survives for more than $7.56$ Gyr.

In the early stages of the galaxy’s evolution, its rotation axis was the minor
axis ($q=\frac{b}{a}<1$). After a massive merger with its satellite halo, the
galaxy’s axis of rotation changes to the major axis. We assumed here
axisymmetry around its major axis ($q=\frac{b}{a}>1$). We created JAM models
for each of the three main stages of the galaxy’s evolution: pre-merger
($z=1.58$), shortly after the merger ($z=0.58$), and at the present day
($z=0.00$).

We used the initial oblate case as a benchmark to evaluate the performance of
our prolate models.

##  3 Methods

The Jeans dynamical models link the generally derived from observations
stellar density and kinematics of the modelled system to its gravitational
potential.

From the simulation, we know the surface brightness profile of the galaxy and
perform a simple Sérsic profile fit to obtain its central surface brightness,
half-light radius, and Sérsic index. We represent the resulting Sérsic profile
with a Multi Gaussian Expansion (MGE), which makes it easy to deproject the
density profile.

We have also created Voronoi maps of the line of sight velocity ($v_{los}$)
and velocity dispersion ($\sigma$) at different snapshots from the simulation.
We assigned uncertainties to these values based on the local surface
brightness (uncertainties increase with decreasing surface brightness).

In our dynamical model, the gravitational potential is driven by DM only. We
used a generalised Navarro-Frank-White (NFW) density profile as representative
of the gravitational potential to solve the Jeans equations. In this
particular JAM model, we used a spherical, cored NFW profile (central density
slope $\gamma=0$), which we later found to give the most promising mass
profile of the model, as compared with the cosmological simulation.

We used Markov Chain Monte Carlo (MCMC) method to get the maximum likelihood
of the fitted parameters. We fit for the central density ($\rho_{0}$) and
scale radius ($r_{0}$) of the NFW profile, as well as velocity anisotropy
$(\beta_{z}=1-\frac{\sigma_{z}^{2}}{\sigma_{r}^{2}})$ and rotation amplitude
($\kappa$) of the galaxy.

##  4 Results

We first modelled the simulated galaxy during its pre-merger stage ($z=1.58$),
when it still had normal oblate rotation. After $2$ Gyr, due to a massive
collision, the galaxy began rotating in a prolate manner. This is the stage in
the simulation, at which we performed our second JAM model ($z=0.58$).
Interestingly, the galaxy also has prolate rotation at $z=0$ (the present day
of the simulation), indicating that it is in a state of quasi-equilibrium.

Figure 1: The top-left figure displays the line-of-sight velocity from the
cosmological simulation of the galaxy at redshift $z=0.00$ in Voronoi bins,
while the top-right figure shows the corresponding JAM model. The black
ellipse indicates the half-light radius. The bottom left figure displays the
velocity dispersion from the cosmological simulation, while the bottom right
shows the dispersion obtained from the fitted JAM model.

In all three cases, we fitted the line-of-sight velocity ($v_{los}$) and
velocity dispersion, and we estimated the mass distribution of the simulated
galaxy using the modelled parameters listed above. An example of the fit for
the last case ($z=0$) is shown in Figure 1. The Voronoi maps of the mock
observations of the galaxy rotation and velocity dispersion are compared to
the best-fit results from the Jeans model.

The Jeans model provided a well matched line-of-sight velocity and velocity
dispersion for all the stages within the half-light radius. Despite the fewer
data points outside the half-light radius, we also got a relatively good fit
with the data.

The capability to determine the galaxy’s mass profile is one of the main
advantages of dynamical modelling. We calculated the mass for each of the
three models stated above. The masses for all three models are displayed in
Figure 2, along with the mass derived from the cosmological simulation.

Figure 2: Best-fit mass profiles of the simulated galaxy at pre-merger
($z=1.58$, left), just after the merger ($z=0.58$, middle), and at the present
time ($z=0$, right), compared to the true mass distribution from the
cosmological simulation. The cross shows the virial mass estimate enclosed
within the half-light radius.

##  5 N-body Simulation of the galaxy

We also checked the survivability of the galaxy’s rotation using an N-body
simulation.

We used the inverse transform sampling method to draw positions from the best-
fit de-projected 3D sérsic profile and initial 3D velocities from the best-fit
Jeans model.

We ran the N-body code assuming star particles move within an embedded
spherical DM halo. The simulation was done using SWIFT (Schaller et al.,
2016). Interestingly, the galaxy loses its prolate rotation after only $100$
Myr, in contrast to its longevity in the cosmological simulation. More
research is necessary to determine the underlying reasons. For example, we
could try changing the shape of the DM halo to prolate as well.

##  6 Conclusions

We successfully applied the JAM code on a cosmologically simulated, prolate
rotating dwarf galaxy and fitted its velocity moments. We used the early
oblate, pre-merger stage as a benchmark for our project. We then modelled the
prolate rotating stages shortly after the merger ($z=0.58$) and at the present
time ($z=0.0$). In all cases, we fitted a spherical NFW mass profile, as well
as the amplitude of rotation and velocity anisotropy. We then calculated the
mass profile of the galaxy at different evolutionary stages and compared them
to the known mass profiles from the simulation, which agree quite well. This
project confirms the viability of the Jeans method for estimating the masses
of prolate rotating galaxies in the approximation of quasi-dynamical
equilibrium. Furthermore, we also did an N-body simulation but in this simple
setup, the galaxy’s prolate rotation dies after $100$ Myr. Further work is
needed to better understand the survivability of its prolate nature.

## Acknowledgements

We want to express our gratitude to Dr. Salvador Cardona-Barrero for providing
the simulation data. We acknowledge Dr. Yves Revaz for his help to run the
N-body simulation. We thank Prof. Dr. Maria-Rosa Cioni for the opportunity to
conduct this research at the AIP.

## References

  * Amorisco et al. (2014) Amorisco, N. C., Evans, N. W., & van de Ven, G. 2014, Nature, 507, 335. doi:10.1038/nature12995
  * Cappellari (2016) Cappellari, M. 2016, Annu. Rev. Astron. Astrophys., 54, 597. doi:10.1146/annurev-astro-082214-122432
  * Cappellari (2020) Cappellari, M. 2020, MNRAS, 494, 4819. doi:10.1093/mnras/staa959
  * Cardona-Barrero et al. (2021) Cardona-Barrero, S., Battaglia, G., Di Cintio, A., et al. 2021, MNRAS, 505, L100. doi:10.1093/mnrasl/slab059
  * Kacharov et al. (2017) Kacharov, N., Battaglia, G., Rejkuba, M., et al. 2017, MNRAS, 466, 2006. doi:10.1093/mnras/stw3188
  * Revaz & Jablonka (2012) Revaz, Y. & Jablonka, P. 2012, A&A, 538, A82. doi:10.1051/0004-6361/201117402
  * Schaller et al. (2016) Schaller, M., Gonnet, P., Chalk, A. B. G., et al. 2016, Proceedings of the Platform for Advanced Scientific Computing Conference, 2. doi:10.1145/2929908.2929916
  * Springel (2005) Springel, V. 2005, MNRAS, 364, 1105. doi:10.1111/j.1365-2966.2005.09655.x


# qecGPT: decoding Quantum Error-correcting Codes with
Generative Pre-trained Transformers

Hanyan Cao  CAS Key Laboratory for Theoretical Physics, Institute of
Theoretical Physics, Chinese Academy of Sciences, Beijing 100190, China
School of Physical Sciences, University of Chinese Academy of Sciences,
Beijing 100049, China  Feng Pan  CAS Key Laboratory for Theoretical Physics,
Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing 100190,
China  Yijia Wang  CAS Key Laboratory for Theoretical Physics, Institute of
Theoretical Physics, Chinese Academy of Sciences, Beijing 100190, China
School of Physical Sciences, University of Chinese Academy of Sciences,
Beijing 100049, China  Pan Zhang<EMAIL_ADDRESS>CAS Key Laboratory for
Theoretical Physics, Institute of Theoretical Physics, Chinese Academy of
Sciences, Beijing 100190, China  School of Fundamental Physics and
Mathematical Sciences, Hangzhou Institute for Advanced Study, UCAS, Hangzhou
310024, China  International Centre for Theoretical Physics Asia-Pacific,
Beijing/Hangzhou, China

###### Abstract

We propose a general framework for decoding quantum error-correcting codes
with generative modeling. The model utilizes autoregressive neural networks,
specifically Transformers, to learn the joint probability of logical operators
and syndromes. This training is in an unsupervised way, without the need for
labeled training data, and is thus referred to as pre-training. After the pre-
training, the model can efficiently compute the likelihood of logical
operators for any given syndrome, using maximum likelihood decoding. It can
directly generate the most-likely logical operators with computational
complexity $\mathcal{O}(2k)$ in the number of logical qubits $k$, which is
significantly better than the conventional maximum likelihood decoding
algorithms that require $\mathcal{O}(4^{k})$ computation.

Based on the pre-trained model, we further propose refinement to achieve more
accurately the likelihood of logical operators for a given syndrome by
directly sampling the stabilizer operators. We perform numerical experiments
on stabilizer codes with small code distances, using both depolarizing error
models and error models with correlated noise. The results show that our
approach provides significantly better decoding accuracy than the minimum
weight perfect matching and belief-propagation-based algorithms. Our framework
is general and can be applied to any error model and quantum codes with
different topologies such as surface codes and quantum LDPC codes.
Furthermore, it leverages the parallelization capabilities of GPUs, enabling
simultaneous decoding of a large number of syndromes. Our approach sheds light
on the efficient and accurate decoding of quantum error-correcting codes using
generative artificial intelligence and modern computational power.

Quantum computers can potentially solve practical problems which are
intractable for classical computers. However, the current implementation of
quantum computers has an issue with noise, which limits its power. An
essential step towards fault-tolerant quantum computing is quantum error
correction (QEC), which now becomes one of the key research frontiers in both
theoretical studies and hardware developments [1, 2] of quantum computation
[3]. In QEC, logical states with $k$ logical qubits are encoded using $n$
physical qubits with redundancy. The effects of continuous errors can be
digitalized into a finite set of discrete errors, which can be obtained by
measuring the redundant ancilla qubits, giving an error syndrome. Then a
decoding algorithm infers the information of errors based on the syndrome and
determines an appropriate operation to correct the logical error. However, the
decoding problem is a hard problem, for example, it belongs to the class of
$\\#$P hard problem in the classical error-correcting codes. In quantum codes,
decoding is considered to be more challenging than classical code, because the
errors inherently degenerate, and the corresponding factor graph for the codes
is more complex, for example, in CSS code the factor graph always contains
loops with various sizes due to the commutation relations, so standard
decoding algorithms such as belief propagation do not work as well as in
classical low-density parity check (LDPC) codes.

While a number of algorithms have been proposed for decoding quantum error-
correcting codes, we lack general decoding algorithms that are efficient and
accurate. The minimum weight perfect matching algorithm [4, 5] can decode
surface code efficiently, however, as a minimum-weight decoder ignores the
degeneracy of quantum codes, in principle its performance usually has a gap to
the theoretical limit. Moreover, it is less efficient in non-planar graphs and
is challenging when applied to code on hypergraphs where the distances between
two nodes are not well defined. As a prototype of the maximum likelihood
decoder (MLD), tensor network methods (e.g. the boundary matrix product state
method [6]) consider the degeneracies of quantum codes and work close to the
theoretical limit in surface code. However, for general codes not defined on
lattices with open boundaries, the tensor network contractions are difficult
to apply due to the large treewidth of the graph. Another issue for existing
maximum-likelihood decoders is computing probabilities for $4^{k}$ logical
operators for $k$ logical qubits, which is intractable for a large $k$.
Moreover, the contraction of tensor networks for each syndrome consumes
significantly more computational resources than the minimum-weight decoding
algorithms and hence less efficient. Recently, a number of neural network
decoders are proposed for leveraging fast inference in neural networks on
modern GPUs [7, 8, 9, 10, 11, 12, 13, 14]. These methods are based on
supervised learning, meaning that training the neural network model requires a
large dataset prepared in advance with labels computed using another teacher
decoding algorithm. Both the dataset size and the accuracy of the teacher
decoding algorithm limit the performance of the supervised neural network
decoders.

In this work, we propose a maximum likelihood decoding approach based on
unsupervised generative modeling in machine learning. The proposed algorithm
enjoys efficient decoding using the fast inference of autoregressive neural
networks, especially on GPUs, the training directly uses the error model and
thus does not require preparing labeled data for training. The autoregressive
neural networks can be applied to quantum codes with arbitrary topology, e.g.
for general quantum low-density parity check (QLDPC) codes, it also supports
directly generating logical operators for an arbitrary number of logical
qubits $k$ by reducing the computational complexity from $\mathcal{O}(4^{k})$
(in computing probabilities of all logical operators in conventional maximum-
likelihood decoding) to $\mathcal{O}(2k)$. In the following text we will first
introduce how to link the decoding of stabilizer codes to the generative
modeling, then introduce the qecGPT, the pre-trained version of our approach
using a specific autoregressive model, casual transformers, then introduce the
refinement of decoding accuracy based on the pre-trained model.

#### Maximum likelihood decoding—

Consider a $[[n,k,d]]$ quantum correction code where a logical state
$\ket{\phi}$ with $k$ logical qubits is encoded using a code word $\ket{\psi}$
with $n$ physical qubits. The minimum distance between the code words is $d$.
When an error occurs on the state $\psi$, it is considered as an effect of
applying an error operator $E$ belonging to the Pauli group
$\mathcal{P}_{n}=\pm i\\{I,X,Y,Z\\}$. In the stabilizer formalism [15, 16],
encoded states are stabilized by some operators $\\{s\\}$, i.e.
$S\ket{\psi}=\ket{\psi}$. The operators form a stabilizer group
$\mathcal{S}=\langle g_{1},g_{n},\cdots,g_{m}\rangle$, which is an Abelian
sub-group of $\mathcal{P}_{n}$, and is generated by $m=n-k$ independent
generators.

When an error $E$ occurs, the encoding state $\psi$ may not be stabilized by
the stabilizers anymore, this can be tested by measuring the ancilla qubits
corresponding to the stabilizer generators, yielding syndrome
$\gamma(E)=\\{\gamma_{1}(E),\gamma_{1}(E),\cdots,\gamma_{m}(E)\\}$, with
$\gamma_{i}(E)=0$ if $g_{i}$ and $E$ commute and $\gamma_{i}(E)=1$ if they
anti-commute. In other words, if the syndrome is not trivial, then the error
$E$ must anti-commute with some of the stabilizer generators. If the syndrome
is trivial, the error $E$ commutates with all $m$ stabilizer generators, then
$E$ is either an element of the stabilizer group or belongs to the logical
operators, which is generated by $\mathcal{L}=\langle
l_{1}^{x},l_{1}^{z},l_{2}^{x},l_{2}^{z},\cdots,l_{k}^{x},l_{k}^{z}\rangle$.
Here $l_{i}^{x}$ and $l_{i}^{z}$ denote the logical X and logical Z operators
of the $i$’th logical qubits respectively. In addition to $2^{m}$ stabilizer
operators and $4^{k}$ logical operators, there are still $2^{m}$ operators
that do not commutate with stabilizer generators, they belong to the pure
error subgroup $\mathcal{E}$ which is Abelian and satisfies the commutation
relation $e_{i}g_{i}=(-1)^{\delta}_{ij}g_{j}e_{i}$. The three subgroups
introduced above indicate a structure of the Pauli group, which is revealed by
the decomposition of the Pauli group
$\mathcal{P}_{n}=\mathcal{E}\otimes\mathcal{L}\otimes\mathcal{S}$.

Based on the $\\{\mathcal{E},\mathcal{L},\mathcal{S}\\}$ decomposition, we can
map an error $E$ to a configuration $\alpha,\beta,\gamma$, and vice versa.
Here $\alpha\in\\{0,1\\}^{m}$ is the configuration of $m$ stabilizer
generators, and each value $\alpha_{i}$ is determined by the commutation
relation between $E$ and the pure error generator $E_{i}$;
$\beta\in\\{0,1\\}^{2k}$ denotes the configuration for logical $X$ and logical
$Z$ operators; $\gamma\in\\{0,1\\}^{m}$ is the configuration of $m$ pure error
generators, and each value $\gamma_{i}$ is determined by the commutation
relation between $E$ and the stabilizer generators. Using this mapping, we can
see the degeneracy of errors, that is, given the logical configuration
$\beta$, there are $2^{m}$ assignments of $\\{\alpha\\}$ give the same
syndrome. So computing likelihood of an logical operator configuration $\beta$
needs to consider all $\alpha$ configurations with

$\displaystyle p(\beta,\gamma)=\sum_{\alpha}p(\alpha,\beta,\gamma).$ (1)

In this sense, we consider the total probability of a coset of the stabilizer
sub-group $\mathcal{S}$, rather than the probability of a single error.

However, there are several challenges for maximum likelihood decoding. The
first one is that computing the coset probability is a $\\#$P problem, no
general exact algorithm exists and approximate algorithms e.g. tensor network
contractions are usually time-consuming; the second challenge is that one
needs to repeat the computation (by summing over all $\alpha$ configurations)
for each syndrome; the third challenge is the exponential computational
complexity in the number of logical qubits $k$ because conventionally, one
needs to enumerate $4^{k}$ logical operators, compute their closet
probabilities, and find the one with the largest probability.

#### Generative maximum likelihood decoding—

We propose to solve the challenges of maximum likelihood decoding using a
framework based on generative modeling. First we approximate the joint
distribution $p(\alpha,\beta,\gamma)$ using a parameterized variational
distribution $q_{\theta}(\alpha,\beta,\gamma)$ satisfying

$q_{\theta}(\alpha,\beta,\gamma)=q(\alpha|\mathbf{\beta},\gamma)q(\beta|\gamma)q(\gamma).$
(2)

Here $q(\alpha|\beta,\gamma)$ is the conditional probability distribution of
stabilizer configuration $\alpha$ given the logical operator $\beta$ and the
syndrome $\gamma$, $q(\beta|\gamma)$ is the conditional probability of $\beta$
given syndrome. $\theta$ denotes the parameters of the variational
distribution $q_{\theta}(\alpha,\beta,\gamma)$. By learning $\theta$, we make
the variational distribution close to the true joint distribution
$p(\alpha,\beta,\gamma)$ given by the error model, and make the variational
conditional distribution $q(\beta|\gamma)$ close to the true conditional
distribution $p(\beta|\gamma)$ which is intractable in general. With an
accurate estimate of conditional probabilities, we can evaluate the likelihood
of logical operators for all syndromes, and generate a configuration of
logical operators by sampling $q(\beta|\gamma)$. In other words, the learned
variational joint distribution satisfies the condition that the configuration
for stabilizer generators $\alpha$ can be traced out automatically.

We further ask the conditional probabilities also satisfy the autoregressive
properties for each variable, i.e.
$q(\beta|\gamma)=q(\beta_{1}|\beta_{2},...,\beta_{2k}|\gamma)q(\beta_{2}|\beta_{3},...,\beta_{2k}|\gamma)\cdots
q(\beta_{2k}|\gamma)$. In this way, all $2k$ logical variables can be
generated one by one following the conditional probabilities. This is known as
ancestral sampling which is an unbiased sampling from the variational
conditional distribution $q(\beta|\gamma)$ [17]. The pictorial representation
of the generative modeling is illustrated in Fig. 1, where we can see that all
the variables are assigned order and each variable only relies on the
variables prior to it. I.e., the conditional probability of configuration of a
variable $s_{i}$ is a function of configurations of variables before it
$\\{s_{1},s_{2},\cdots,s_{i-1}\\}=\mathbf{s}_{<i}$, with
$q(s_{i}|\mathbf{s}_{<i})$. This property of the parameterization is known as
the autoregressive property, also known as causal property, if we regard the
variables before it as its “history”, and the variables behind it as its
“future”. Many neural network models satisfy this property and are known as
autoregressive neural networks, especially in the models for natural languages
where the words are generated one by one.

#### The Generative Pre-trained Transformers—

In this work, to parameterize $q_{\theta}(\alpha,\beta,\gamma)$ we adopt the
Transformers, one of the most powerful autoregressive neural networks [18] and
has been used in many applications including chatGPT [19, 20, 21]. We use the
decoder layer of the Transformer composed of an embedding layer and a
positioning encoding layer which map the input configuration to a higher
dimensional feature space, attention layers with triangular mask (to ensure
the autoregressive properties), and a linear layer and a sigmoid output layer
that outputs the joint probability distributions. The details of the
Transformer can be found in the appendix.

Figure 1: Illustration of the structure of qecGPT.

Figure 2: Logical error rates of our algorithm (qecGPT and qecGPT+refinement)
with different physical error rates, compared with MWPM, and BPOSD algorithm
on the [13, 1, 3] surface code (left), [41, 1, 5] surface code (middle), and
[12, 1, 2] 3D-surface code (right). The error model is the depolarizing model.
Each data point in the figures is averaged over $10000$ error instances. The
black lines are the optimal maximum-likelihood decoding algorithm which
exactly sums all stabilizer configurations. The insets are the difference
between the approximate algorithms and the exact algorithm.

The parameters are learned to minimize the distance between the true
distribution $p(\alpha,\beta,\gamma)$ and the variational distribution
parameterized using the Transformer. In this work, we assume that we have
samples of the noise model $p(E)$, we choose the forward Kullback-Leibler
divergence as the distance measure of two probability distributions

$\displaystyle
D_{\mathrm{K}L}(p|q)=\sum_{\alpha,\beta,\gamma}p(\alpha,\beta,\gamma)\log\frac{p(\alpha,\beta,\gamma)}{q_{\theta}(\alpha,\beta,\gamma)}.$

This yields a negative log-likelihood loss function

$\hat{\theta}=\mathop{\arg\min}_{\theta}D_{\mathrm{K}L}=\mathop{\arg\min}_{\theta}\left[-\sum_{\alpha,\beta,\gamma\sim
p}\log q_{\theta}(\alpha,\beta,\gamma)\right].$ (3)

And the parameters are updated using a gradient-based optimizer. After
training, given a syndrome $\gamma$, we can generate a configuration of the
logical operator $\beta=\\{\beta_{1},\beta_{2},\cdots,\beta_{2k}\\}$ one by
one using the learned conditional probabilities

$\hat{\beta}_{i}=\mathop{\arg\max}_{\beta_{i}}q(\beta_{i}|\beta_{1},\beta_{2},\cdots,\beta_{i-1},\gamma_{1},\gamma_{2},\cdots,\gamma_{m}).$

So a logical configuration is generated variable-by-variable given a syndrome.
This is analogous to the generation of text from the chatGPT [19, 20, 21],
where the text is generated word-by-word given a prompt. Notice that by
training once, a single variational distribution gives the conditional
probabilities of all $2^{m}$ syndromes, maximizing the likelihood of all
syndromes. So we term it pre-training. The advantage of the pre-training is
that the conditional probability for any syndrome can be computed efficiently
using a single pass of the neural network. We can further optimize the
accuracy of the conditional probability based on the pre-trained model given a
particular syndrome, which we term as refinement. There could be approaches
for the refinement, for example, we can minimize the
$D_{\mathrm{K}L}(q_{\theta}(\alpha,\beta|\gamma)|p(\alpha,\beta|\gamma))$
using e.g. the method of variational autoregressive networks. In this work, we
propose a straightforward way for the refinement with a small number of
logical qubits $k$, taking advantage of the generative modeling. In addition
to generating $\beta$ configurations, we can also generate efficiently the
stabilizer configurations $\alpha$ using $q(\alpha,\beta,\gamma)$, and use
them to evaluate an unbiased estimate of the joint probability of $\beta$ and
$\gamma$

$\displaystyle p(\beta,\gamma)$
$\displaystyle=\sum_{\alpha}q_{\theta}(\alpha,\beta,\gamma)\frac{p(\alpha,\beta,\gamma)}{q_{\theta}(\alpha,\beta,\gamma)}\approx\frac{1}{N}\sum_{\alpha\sim
q}\frac{p(\alpha,\beta,\gamma)}{q_{\theta}(\alpha,\beta,\gamma)}.$

Here we use the samples of the variational distribution and the reweighting to
compute an unbiased estimate of the joint distribution, $N$ is the number of
samples.

An advantage of our approach is the insensitivity to the topology of the code,
i.e. the connectivity of the stabilizer generators, or in other words, the
structure of the parity check matrix. The transformer representation of the
variational distribution can be used for any code topology without modifying
the structure of the transformer, thanks to the self-attention mechanism which
can automatically capture correlations in variables. Our approach is also
insensitive to the true parameters used in the error model. For example, under
the depolarizing noise model, we can train qecGPT with a particular physical
error rate smaller than the threshold, and use our model to decode error under
the depolarizing noise model with distinct error rates, without significantly
increasing the logical error rate. This effect is quite common in inference
problems with mismatched parameters, e.g. in [22]. We refer to the Appendices
for details.

#### Numerical experiments—

We evaluate our algorithm by comparing the logical error rate of our algorithm
to the minimum weight perfect matching (MWPM) [5] and belief propagation
augmented by ordered statistics decoding (BPOSD) [23] algorithm on the surface
code. In Fig. 2 the surface codes have $k=1$ logical qubits. In Fig. 2(left)
the distance $d=3$ is small, and we see that the pre-trained model qecGPT
performs very close to the exact maximum likelihood decoding which computes
the exact likelihood for all $4^{k}$ logical operators using exact tensor
network contractions. In Fig. 2 middle, the system is larger with $d=5$, and
we see that the pre-trained model gives slightly worse results than the exact
algorithm while still much better than MWPM and BPOSD. We also see that the
refinement significantly improves the performance of qecGPT, making the
performance very close to the optimal MLD decoder.

To demonstrate the generality of our approach with different code topologies,
we also test the algorithm in the stabilizer code on a 3-dimensional lattice,
which is usually termed as 3D surface code [24]. Note that for the 3D
structure code, the MWPM algorithm does not apply directly, so we only compare
the logical error rates of qecGPT with BPOSD. We can see from Fig. 2 that the
logical error of qecGPT coincides very well with the exact MLD algorithm and
significantly outperforms BPOSD. Additional numerical results and performance
comparisons with $k>1$ logical qubits and under a noise model with correlated
noise can be found in Appendices.

#### Discussions—

We have introduced a general framework for decoding quantum error correction
code with generative modeling. Our method approximates the joint distribution
of errors using variational autoregressive neural networks. We propose a pre-
trained model for the fast generation of the maximum-likelihood logical
operators and a refinement to increase the accuracy given a syndrome. The
advantage of our generative molding is that it solves the difficulties of the
maximum likelihood decoding in summing over $2^{n}$ stabilizer configurations
and in computing the probability of $4^{k}$ logical configurations. Another
advantage is its generality in code topologies, e.g. it can be applied to 2D
codes and QLDPCs code without modifying the model or the algorithm.

In our work, we have successfully trained the qecGPT using a single GPU and
conducted experiments on small codes with distances up to $d=7$. Although the
decoding process is fast and efficient, the training phase is slow and poses a
challenge when it comes to applying it to a larger code. However, we believe
that this bottleneck can be resolved by using a larger model and exploiting
more computational resources, such as multiple GPUs or even a supercomputer.
This approach is similar to the behavior of chatGPT, which has shown
remarkable performance when trained with a large amount of data and
computational power, as reported in [21, 19]. We intend to explore this avenue
in the future and see how it can further improve the qecGPT’s efficiency and
scalability.

###### Acknowledgements.

A python implementation and a Jupyter Notebook tutorial of our algorithm are
available at [25]. We thank Weilei Zeng, Lingling Lao, and Ying Li for their
helpful discussions and Michael Vasmer for providing 3D surface code data.

## References

  * Panteleev and Kalachev [2022a] P. Panteleev and G. Kalachev, Quantum LDPC codes with almost linear minimum distance, IEEE Transactions on Information Theory 68, 213 (2022a).
  * Panteleev and Kalachev [2022b] P. Panteleev and G. Kalachev, Asymptotically good quantum and locally testable classical ldpc codes, in  _Proceedings of the 54th Annual ACM SIGACT Symposium on Theory of Computing_, STOC 2022 (Association for Computing Machinery, New York, NY, USA, 2022) p. 375–388.
  * Google Quantum AI [2023] Google Quantum AI, Suppressing quantum errors by scaling a surface code logical qubit, Nature 614, 676 (2023).
  * Dennis  _et al._ [2002] E. Dennis, A. Kitaev, A. Landahl, and J. Preskill,  _Topological quantum memory_, Tech. Rep. 9 (2002).
  * Higgott [2021] O. Higgott, Pymatching: A python package for decoding quantum codes with minimum-weight perfect matching (2021), arXiv:2105.13082 [quant-ph] .
  * Bravyi  _et al._ [2014] S. Bravyi, M. Suchara, and A. Vargo, Efficient algorithms for maximum likelihood decoding in the surface code, Physical Review A 90, 032326 (2014).
  * Torlai and Melko [2017] G. Torlai and R. G. Melko, Neural decoder for topological codes, Phys. Rev. Lett. 119, 030501 (2017).
  * Varsamopoulos  _et al._ [2017] S. Varsamopoulos, B. Criger, and K. Bertels, Decoding small surface codes with feedforward neural networks, Quantum Science and Technology 3 (2017).
  * Krastanov and Jiang [2017] S. Krastanov and L. Jiang, Deep Neural Network Probabilistic Decoder for Stabilizer Codes, Scientific Reports 7 (2017).
  * Varsamopoulos  _et al._ [2020] S. Varsamopoulos, K. Bertels, and C. G. Almudever, Comparing neural network based decoders for the surface code, IEEE Transactions on Computers 69, 300 (2020).
  * Overwater  _et al._ [2022] R. W. Overwater, M. Babaie, and F. Sebastiano, Neural-Network Decoders for Quantum Error Correction Using Surface Codes: A Space Exploration of the Hardware Cost-Performance Tradeoffs, IEEE Transactions on Quantum Engineering 3, 1 (2022).
  * Baireuther  _et al._ [2018] P. Baireuther, T. E. O’Brien, B. Tarasinski, and C. W. Beenakker, Machine-learning-assisted correction of correlated qubit errors in a topological code, Quantum 2, 48 (2018).
  * Davaasuren  _et al._ [2020] A. Davaasuren, Y. Suzuki, K. Fujii, and M. Koashi, General framework for constructing fast and near-optimal machine-learning-based decoder of the topological stabilizer codes, Physical Review Research 2 (2020).
  * Gicev  _et al._ [2021] S. Gicev, L. C. Hollenberg, and M. Usman, A scalable and fast artificial neural network syndrome decoder for surface codes, arXiv preprint arXiv:2110.05854 (2021).
  * Gottesman [1997] D. Gottesman,  _Stabilizer codes and quantum error correction_ (California Institute of Technology, 1997).
  * Nielsen and Chuang [2010] M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information: 10th Anniversary Edition_ (Cambridge University Press, 2010).
  * Bishop [2006] C. M. Bishop,  _Pattern Recognition and Machine Learning (Information Science and Statistics)_ (Springer-Verlag, Berlin, Heidelberg, 2006).
  * Vaswani  _et al._ [2017] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin, Attention is all you need, Advances in neural information processing systems 30 (2017).
  * Radford  _et al._ [2018] A. Radford, K. Narasimhan, T. Salimans, and I. Sutskever, Improving language understanding by generative pre-training (OpenAI, 2018).
  * [20] https://openai.com/chatgpt.
  * OpenAI [2023] OpenAI, Gpt-4 technical report (2023), arXiv:2303.08774 [cs.CL] .
  * Zhang and Moore [2014] P. Zhang and C. Moore, Scalable detection of statistically significant communities and hierarchies, using message passing for modularity, Proceedings of the National Academy of Sciences of the United States of America 111, 18144 (2014).
  * Roffe  _et al._ [2020] J. Roffe, D. R. White, S. Burton, and E. Campbell, Decoding across the quantum low-density parity-check code landscape, Physical Review Research 2, 10.1103/PhysRevResearch.2.043423 (2020).
  * Vasmer and Browne [2019] M. Vasmer and D. E. Browne, Three-dimensional surface codes: Transversal gates and fault-tolerant architectures, Phys. Rev. A 100, 012312 (2019).
  * [25] https:github.com/CHY-i/qecGPT.
  * Kolmogorov [2009] V. Kolmogorov, Blossom V: A new implementation of a minimum cost perfect matching algorithm, Mathematical Programming Computation 1, 43 (2009).
  * Fakoor  _et al._ [2020] R. Fakoor, P. Chaudhari, J. Mueller, and A. J. Smola, Trade: Transformers for density estimation (2020), arXiv:2004.02441 [cs.LG] .
  * Williams [1992] R. J. Williams, Simple statistical gradient-following algorithms for connectionist reinforcement learning, Machine learning 8, 229 (1992).
  * Wu  _et al._ [2019] D. Wu, L. Wang, and P. Zhang, Solving statistical mechanics using variational autoregressive networks, Phys. Rev. Lett. 122, 080602 (2019).
  * Bombin  _et al._ [2012] H. Bombin, R. S. Andrist, M. Ohzeki, H. G. Katzgraber, and M. A. Martin-Delgado, Strong resilience of topological codes to depolarization, Physical Review X 2, 10.1103/physrevx.2.021004 (2012).
  * Battistel  _et al._ [2023] F. Battistel, C. Chamberland, K. Johar, R. W. Overwater, F. Sebastiano, L. Skoric, Y. Ueno, and M. Usman, Real-time decoding for fault-tolerant quantum computing: Progress, challenges and outlook, arXiv preprint arXiv:2303.00054 (2023).
  * Bombin and Martin-Delgado [2007] H. Bombin and M. A. Martin-Delgado, Optimal resources for topological two-dimensional stabilizer codes: Comparative study, Physical Review A 76 (2007).

##  Appendix A Stabilizer codes

The stabilizer code [15] is a very important class of quantum error-correcting
codes. Here we will first describe the stabilizer codes and then introduce the
decoding algorithms. Consider a $[n,k,d]$ stabilizer code, states of $k$
logical qubits are encoded to $n$ physical qubits states. The states
$\ket{\phi}$ of $n$ physical form a $2^{n}$ Hilbert space $\mathcal{H}_{n}$,
thus the encoding states $\ket{\psi}$ of $k$ logical qubits form a subspace of
$\mathcal{H}_{n}$ and can be represented by superposition of $\ket{\phi}$. The
bit-flip $X$ and phase-flip $Z$ errors may occur on a single qubit state.
Then, for $n$ qubits state $\ket{\psi}\in\mathcal{H}_{n}$, all errors form a
group $\mathcal{P}_{n}=\mathcal{P}^{\otimes n}$ Pauli group. An elements
$E\in\mathcal{P}_{n}$ acting on $\ket{\psi}$ may cause an error state
$\ket{\psi^{\prime}}$. The quantum error correction is to find a recover
operator $E^{\prime}$ to correct the state
$E^{\prime}\ket{\psi^{\prime}}=\ket{\psi}$. A straightforward idea is to find
which error $E$ has occurred and $E^{\prime}=E$ because of the self-inverse
property of Pauli operators. However, a special encoding allows us to find
only a collection of operators and any operator belonging to this collection
can recover the error. Such an encoding method is called stabilizer code due
to the construction is based on a subgroup $\mathcal{S}$ of $\mathcal{P}_{n}$
called stabilizer group. This group satisfies the following properties:

(a). $\mathcal{S}$ is an abelian group.

(b). $-I\notin\mathcal{S}$.

Then the encoding states can be chosen as follow:

$\centering\\{\ket{\psi}\ \mid\ S\ket{\psi}=\ket{\psi},\ \forall
S\in\mathcal{S},\ \forall\ket{\psi}\in\mathcal{H}_{n}\\}\@add@centering$ (4)

According to the properties of $\mathcal{S}$, there are $2^{m}$ elements in
$\mathcal{S}$, where $m$ denotes number of generators $\langle
g_{1},\cdots,g_{m}\rangle$. If there exists a error $E\in\mathcal{P}_{n}$, one
may observe $\gamma(E)$ with length m called error syndrome

$\centering\gamma(E)_{i}=\left\\{\begin{split}&0,\quad[g_{i},E]=0\\\
&1,\quad\\{g_{i},E\\}=0\\\ \end{split}\right.\@add@centering$ (5)

The stabilizers are usually described by a parity check matrix $H$ with size
$m\times 2n$, where $m=n-k$ denotes the number of stabilizer generators.

####  $\mathbf{F_{2}}$ representaiton—

The $\mathbf{F}_{2}$ representation is an isomorphic map from the Pauli group
to itself. Under the $\mathbf{F}_{2}$ representation, single qubit Pauli
operators are represented by two binary numbers:

$I\rightarrow 00\quad X\rightarrow 01\quad Z\rightarrow 10\quad Y\rightarrow
11$ (6)

At this time, any n-qubits Pauli operator is represented by a binary vector
with length $2n$ [15]. The group multiplication is just the addition(mod 2)
between vectors. And the commutation relation between two operators $A$ and
$B$ can be represented by:

$\centering\mathbf{A}\cdot\Lambda\cdot\mathbf{B}^{T}=\left\\{\begin{split}&0,\quad[A,B]=0\\\
&1,\quad\\{A,B\\}=0\\\ \end{split}\right.\@add@centering$ (7)

In Eq. 7 the bold letters indicate $\mathbf{F}_{2}$ representations of
operators. The dot symbol $\cdot$ denotes matrix multiplication. And the
$\Lambda$ is a $2n\times 2n$ matrix $\begin{pmatrix}0&I\\\ I&0\\\
\end{pmatrix}$. Further, the generators of the stabilizer group form a
$m\times 2n$ matrix $H$ named parity check matrix.

We have another two subgroups of $\mathcal{P}_{n}$. One is the pure errors
group $\mathcal{E}$. The group $\mathcal{E}$ is an Abelian group with $2^{m}$
elements and all generators $\langle e_{1},\cdots,e_{m}\rangle$ satisfy

$\left\\{\begin{split}e_{i}g_{j}&=(-1)^{\delta_{ij}}g_{j}e_{i}\\\
e_{i}e_{j}&=e_{j}e_{i}\end{split}\right.$ (8)

All generators of the pure error group can be stored as another matrix $M_{E}$
with size $m\times 2n$, satisfying

$\begin{split}&H\cdot M_{E}^{T}=I_{m\times m}\\\ &H\cdot
M_{E}^{L}=0.\end{split}$ (9)

And it can be determined with a matrix $D^{\prime}$ which is further computed
using Gaussian elimination on matrix $D=(H|I_{m\times m})$. $D^{\prime}$ can
be organized as $D^{\prime}=(A|B)$, where matrix $A$ is a row echelon matrix.
And each row of $M_{E}$ can be solved by these new equations:

$A\cdot(M_{E})_{i}^{T}=B_{i}$ (10)

Since the number of rows of matrix $H$ is less than the number of columns,
there are some free variables. For simplicity, we fix these variables to $0$.
Moreover, we want these pure error generators to commute to each other as
defined in Eq. (8).

$M_{E}\cdot\Lambda\cdot M_{E}^{T}=0$ (11)

This requires that some stabilizer generators acting on these operators are
determined. It is equivalent to adding some rows of $H$, corresponding to the
stabilizer generators, on rows of $M_{E}$. The pseudo-code for the whole
process is given as follows: [H] Find Pure Errors $M_{E}$

0: Parity Check Matrix $H$

0: $M_{E}$

$m\times 2n=|H|$

$D=(H|I_{m\times m})$

$D^{\prime}(A|B)=GE(D)$, where A is a row echelon matrix.

$M_{E}\leftarrow$ Each row of $M_{E}$ can be solved from the equation
$\sum_{j}A_{kj}(M_{E}^{T})_{ji}=B_{ki}$. All free variables are set to 0.

for $i\in[1,m]$ do

for $j>i$ do

$s_{ij}=(M_{E})_{i}\cdot\Lambda\cdot(M_{E})_{j}^{T}$

if $s_{ij}\neq 0$ then

$(M_{E})_{i}=(M_{E})_{i}+H_{j}$

else

$(M_{E})_{i}=(M_{E})_{i}$

end if

end for

end for

return $M_{E}$

In addition to the stabilizer group and the pure error group, another sub-
group is the logical-operator group $\mathcal{L}$. It represents the logical
errors of logical qubits, is a non-Abelian group, and is generated by $\langle
l_{x1},l_{z1}\cdots,l_{xk},l_{zk}\rangle$ satisfying

$\left\\{\begin{split}&l_{(x/z)i}g_{j}=g_{j}l_{(x/z)i}\\\
&l_{(x/z)i}e_{j}=e_{j}l_{(x/z)i}\\\
&l_{xi}l_{zj}=(-1)^{\delta_{ij}}l_{zj}l_{xi}\end{split}\right.$ (12)

This group has $4\times 2^{2k}$ elements, the constant $4$ comes from the
overall phase $\\{\pm{1},\pm{i}\\}$. However, during the actual error
correction process the overall phase is always ignored. The matrix of
generators of the logical subgroup, $M_{L}$, can be determined given $H$ and
$M_{E}$. As defined in Eq. 12, i.e., the $M_{L}$ is the kernel of matrix
$M=\left(\begin{array}[]{cc}&H\\\ &M_{E}\end{array}\right)$. Here $M$ is a
matrix with size $2m\times 2n$, Gaussian elimination of $M$ gives $M^{\prime}$
and there are $2k$ free variables. The $M_{L}$ is actually a set of bases of
the kernel space, each two rows of the $M_{L}$ must be linearly independent.
Therefore we choose these free variables as a one-hot vector
$(0,\cdots,1_{i},\cdots,0)$ for $i$th row of $M_{L}$. And for satisfying the
condition Eq. 12. We can perform the symplectic Gram-Schmidt orthogonalization
procedure (SGSOP) on $M_{L}$, finding $k$ pairs $L_{x}$ and $L_{z}$. The
pseudo-code is described in Algorithm A. [H] Find Logical Operators $M_{L}$

0: $H$, $M_{E}$

0: $M_{L}$

$M=\left(\begin{array}[]{cc}&H\\\ &M_{E}\end{array}\right)$

$|M|=2m\times 2n$

$M^{\prime}=GE(M)$, where $M^{\prime}$ is a row echelon matrix.

$M_{L}\leftarrow$ Each row of $M_{E}$ can be solved from the equation
$M^{\prime}\cdot(M_{L}^{T})_{i}=0$. The free variables are set to
$(0_{1},\cdots,1_{i},\cdots,0_{2k})$ for $i$th row of $M_{L}$.

for $i\in[1,k]$ do

for $j\in[i+1,2k]$ do

if $(M_{L})_{i}\cdot\Lambda\cdot(M_{L})_{j}^{T}\neq 0$ then

$(M_{L})_{i+1}\leftrightarrow(M_{L})_{j}$

break

end if

end for

end for

for $i\in[1,k]$ do

for $j\in[i+2,2k]$ do

if $(M_{L})_{i}\cdot\Lambda\cdot(M_{L})_{j}^{T}\neq 0$ then

$(M_{L})_{j}=(M_{L})_{j}+(M_{L})_{i+1}$

else if $(M_{L})_{i+1}\cdot\Lambda\cdot(M_{L})_{j}^{T}\neq 0$ then

$(M_{L})_{j}=(M_{L})_{j}+(M_{L})_{i}$

end if

end for

end for

return $M_{L}$

This is a general algorithm to find logical operators for a given $H$.
Actually, one can not distinguish the logical X and logical Z through this
algorithm. However for the CSS code, the parity check matrix of which can
always be written as $H=(H_{z}|H_{x})$, the $H_{x}$ and $H_{z}$ are always
treated respectively.

##  Appendix B The $\\{\mathcal{E}LS\\}$ decomposition

Note that under the $\mathbf{F_{2}}$ representation, any two operators are
commutative to each other. And the information of anticommutation is stored in
a special class of inner products Eq. 7. This allows us to use a more
efficient way to represent Pauli operators. The Pauli group becomes a self-
inverse and Abelian group, under the $\mathbf{F}_{2}$ representation. Thus any
error operator can be generated by generators $\langle
e_{1},\cdots,e_{m},l_{x1},l_{z1}\cdots,l_{xk},l_{zk},g_{1},\cdots,g_{m}\rangle$
and their powers $(\gamma,\beta,\alpha)$.

$E=\prod_{i,j,k}e^{\gamma_{i}}_{i}\times l^{\beta_{j}}_{j}\times
g^{\alpha_{k}}_{k},$ (13)

with
$\alpha=\\{\alpha_{k}\\}\in\\{0,1\\}^{m},\,\,\,\beta=\\{\beta_{j}\\}\in\\{0,1\\}^{2k},\,\,\,\gamma=\\{\gamma_{i}\\}\in\\{0,1\\}^{m},\,\,\,E\in\mathcal{P}_{n}$.
This means that there is a correspondence between an error operator $E$ and a
$\alpha,\beta,\gamma$ configuration.

$E\Longleftrightarrow(\alpha,\beta,\gamma)$ (14)

We term the power-configuration $(\gamma,\beta,\alpha)$ the $\mathcal{E}LS$
configuration of an operator which forms a binary vector with length $2n$
$\\{\alpha,\beta,\gamma\\}^{2n}$. Given an $\mathcal{E}LS$ configuration, one
only needs a series of vector additions (under the $\mathbf{F}_{2}$
representation) to generate the corresponding operator. On the other hand,
given an error operator, the corresponding $\mathcal{E}LS$ configuration can
be determined using Eq. 7 and

$\centering\alpha(E)_{i}=\left\\{\begin{split}&0,\quad[e_{i},E]=0\\\
&1,\quad\\{e_{i},E\\}=0\\\ \end{split}\right.\@add@centering$ (15)

$\centering\beta(E)_{(x/z)i}=\left\\{\begin{split}&0,\quad[l_{(z/x)i},E]=0\\\
&1,\quad\\{l_{(z/x)i},E\\}=0\\\ \end{split}\right.\@add@centering$ (16)

$\centering\gamma(E)_{i}=\left\\{\begin{split}&0,\quad[g_{i},E]=0\\\
&1,\quad\\{g_{i},E\\}=0\\\ \end{split}\right.\@add@centering$ (17)

##  Appendix C The minimum weight decoder

Decoding is to determine the recovery operator of the quantum error correction
code given the syndrome. There are basically two kinds of decoding algorithms.
The first kind of decoding algorithm is known as the minimum weight decoder,
it determines an error with the maximum probability that satisfies all the
syndrome constraints. In this sense, when a nontrivial syndrome $\gamma$ has
been measured, the minimum weight decoding algorithm finds an error operator
$\hat{E}(\gamma)=\mathop{\arg\max}_{E}P(E(\gamma))$, where $P(E(\gamma)))$ is
the generation probability in the error model which is consistent with the
syndrome. The most famous algorithm of the first kind is the Minimum Weight
Perfect Matching (MWPM) algorithm, which assigns a weight to each edge in the
code graph using the probability of the error in the error model [5], then
finds the shortest error chain given the syndrome, it can be done by employing
an efficient algorithm e.g. Blossom algorithm [26]. The MWPM algorithm can
decode the surface code efficiently. However, it is challenging for the MWPM
algorithm to decode when the code graph is a hypergraph, where each edge of
the code graph links to more than two nodes and the distance between two nodes
is ill-defined [5]. The main limitation of the minimum weight decoding
algorithm is that the error with the maximum probability may not be the right
recovery operator because of the degeneracy of the quantum code.

##  Appendix D The maximum likelihood decoder

The second kind of decoding algorithm is known as the maximum likelihood
decoder. Notice that in the $\mathcal{E}LS$ decomposition, the stabilizers,
and the logical operators do not modify the syndrome $\gamma$. This means that
any element of the normalizer of the stabilizer group,
$\mathcal{N}(S)=\mathcal{L}\otimes\mathcal{S}$, does not change the syndrome.
One can imagine that given an error that produces the syndrome, one can apply,
on top of the error, any element from $\mathcal{N}(s)$ without modifying the
syndrome. So in principle, instead of considering a single error that produces
the syndrome, one should consider all possible errors that produce the
syndrome, which form a closet of the $\mathcal{N}(S)$. The maximum likelihood
decoder determines a logical operator $L$ by considering all operators in the
same equivalent class $C(L,\gamma)$ rather than considering a single error,
equivalently, it is summing the probabilities of the closet
$C(\beta(L),\gamma)$, where $\beta(L)$ is the $\beta$ configuration
corresponding to the logical operator $L$. This respects the degeneracy of
quantum codes and is the best decoding algorithm one could do. Obviously, the
coset probability can be computed by summing probabilities of all elements of
the stabilizer group with a particular syndrome and logical operator. The
summation can be done by considering all possible $\alpha=\\{0,1\\}$
configurations.

$\displaystyle\hat{L}=\mathop{\arg\max}_{L\in\mathcal{L}}P(C(L,\gamma))=\sum_{S\in\mathcal{S}}P(E(\gamma)\times
L\times S)=\sum_{\mathbf{\alpha}}P(e(\gamma)\times
L\times(g^{\alpha_{1}}_{1}\times\cdots\times g^{\alpha_{m}}_{m})),$ (18)

where $\times$ denotes the multiplication of group elements, $e(\gamma)$ is
the pure error corresponding to syndrome the Using the representation of
$\beta$ configuration for a logical operator, we have

$\displaystyle\hat{\beta}$
$\displaystyle=\mathop{\arg\max}_{\beta}\sum_{\alpha}P(\alpha,\beta,\gamma).$
(19)

The computation of summing all possible $\alpha$ configurations is analogous
to the computation for the partition function of an Ising spin glass, but one
needs to do the computation for all possible $\beta$ configurations. This
computation belongs to the #P problem and there is no exact algorithm to solve
it in general in polynomial time. The exact computation of the maximum likely
logical operator requires an exponential algorithm. For some special cases,
e.g. code in the 2D lattice such as the Surface code, the summation can be
approximately computed using the tensor network contractions (e.g. with the
boundary matrix product states method) [6] of a 2D tensor network constructed
for a given syndrome, which could be time-consuming. In addition to the
computational cost, the tensor network contraction method has several
limitations. The first limitation is that it is difficult to generalize to
code on other topologies such as on 3-dimensional lattice, or codes with long-
range interactions as in qLDPC, due to the fast increase of the computational
cost. and the decrease in accuracy with a topology having large treewidth; the
second limitation is that one needs to perform tensor network contractions for
each $\beta$ configuration.

##  Appendix E The generative decoder with Transformers

As described in the main text, our generative decoder models the joint
distribution $P(\alpha,\beta,\gamma)$ using the autoregressive neural
networks, it factorizes the joint distribution as a product of conditional
distributions

$q_{\theta}(\alpha,\beta,\gamma)=q(\alpha|\beta,\gamma)q(\beta|\gamma)q(\gamma)$
(20)

Note that changing the order of the $\alpha,\beta,\gamma$ variables does not
change the joint probability distribution due to the Abelian nature of the
subgroups. In this work, we always put the logical variable $\beta$ in the
middle as in Eq. 20. The benefit of this is that the $\beta$ can be always
determined prior to the stabilizer variables $\alpha$ in the joint
distribution. In this way, decoding is performed by sampling the $\beta$
variables from the marginal distribution

$q(\beta,\gamma)=\sum_{\alpha}q_{\theta}(\alpha,\beta,\gamma).$

Furthermore, we ask the conditional probabilities are organized in such as way
$p(\beta_{i}|\beta_{j<i})$ that allows generating a $\beta$ configuration
among all $4^{k}$ possible ones variable-by-variable since we have already
stored all conditional probabilities for each $\beta_{i}$ variables. In this
way, we can reduce the computational complexity of generating a maximum-
likelihood logical operator out of $4^{k}$ logical operators to
$\mathcal{O}(2k)$.

Figure 3: A pictorial illustration of decoding of the quantum error correction
code with a generative pretrained transformer.

The autoregressive model we used here for representing
$q_{\theta}(\alpha,\beta,\gamma)$ is the encoder part of a Transformer [18,
27] with a mask to ensure the autoregressive property. It is also known as
causal transformer, the structure is shown in Fig. 1. The input of the
Transformer is the configuration $(\alpha,\beta,\gamma)$. The embedding layer
increases the dimension of input to the dimension of a model. The information
on position will be added and learned from the position encoding. A triangular
mask is added in the attention block before the Softmax layer to ensure that
each conditional probability of variable $i$ only depends on the variables
before $i$ in the input configuration. Multiple transformer encoder layers are
added after position encoding, and the final Linear layer maps the data from
the model dimension to a length $n$ vector which is the same as the input. The
output is a vector $(\hat{\alpha},\hat{\beta},\hat{\gamma})$, which uses
Sigmoid functions to represent the Bernoulli distributions for the conditional
probabilities. For example, as illustrated in Fig. 1, we have

$\displaystyle\hat{\gamma}_{1}$
$\displaystyle=\sigma(F_{\gamma_{1}}(\gamma_{1}))=q(\gamma_{1})$
$\displaystyle\hat{\gamma}_{2}$
$\displaystyle=\sigma(F_{\gamma_{2}}(\gamma_{1}))=q(\gamma_{2}|\gamma_{1})$
$\displaystyle\hat{\beta}_{1}$
$\displaystyle=\sigma(F_{\beta_{1}}(\gamma_{1},\gamma_{2}))=q(\beta_{1}|\gamma_{1},\gamma_{2})$
$\displaystyle\hat{\beta}_{2}$
$\displaystyle=\sigma(F_{\beta_{2}}(\beta_{1},\gamma_{1},\gamma_{2}))=q(\beta_{2}|\beta_{1},\gamma_{1},\gamma_{2})$
$\displaystyle\hat{\alpha}_{1}$
$\displaystyle=\sigma(F_{\alpha_{1}}(\alpha_{1},\beta_{1},\gamma_{1},\gamma_{2}))=q(\alpha_{1}|\beta_{2},\beta_{1},\gamma_{1},\gamma_{2})$
$\displaystyle\hat{\alpha}_{2}$
$\displaystyle=\sigma(F_{\alpha_{2}}(\alpha_{2},\alpha_{1},\beta_{1},\gamma_{1},\gamma_{2}))=q(\alpha_{2}|\alpha_{1},\beta_{2},\beta_{1},\gamma_{1},\gamma_{2})$
(21)

where $\sigma(\cdot)$ is the Sigmoid function, and function $F$ denotes the
map of the neural network. We can see that the product of all the output of
the transformer gives the joint distribution, as

$\displaystyle\hat{\alpha}_{2}\hat{\alpha}_{1}\hat{\beta}_{2}\hat{\beta}_{1}\hat{\gamma}_{2}\hat{\gamma}_{1}$
$\displaystyle=q(\alpha_{2}|\alpha_{1},\beta_{2},\beta_{1},\gamma_{2},\gamma_{1})q(\alpha_{1}|\beta_{2},\beta_{1},\gamma_{2}\gamma_{1})q(\beta_{2}|\beta_{1},\gamma_{2},\gamma_{1})q(\beta_{1}|\gamma_{2},\gamma_{1})q(\gamma_{2}|\gamma_{1})q(\gamma_{1})$
$\displaystyle=q(\alpha_{2},\alpha_{1},\beta_{2},\beta_{1},\gamma_{2},\gamma_{1}).$
(22)

Another important property (but not obvious) that we can obtain from the
product of the output of the Transformer is the marginal distribution.

$\displaystyle\hat{\beta}_{2}\hat{\beta}_{1}$
$\displaystyle=q(\beta_{2}|\beta_{1},\gamma_{2},\gamma_{1})q(\beta_{1}|\gamma_{2},\gamma_{1})$
$\displaystyle=q(\beta_{2},\beta_{1}|\gamma_{2},\gamma_{1}).$ (23)

Which gives the normalized conditional probability of logical variables given
a syndrome, and surprisingly ignores the stabilizer configurations $\alpha$.
Based on this conditional probability one can evaluate the likelihood of the
logical operators and perform the decoding after the neural network is well
trained.

##  Appendix F Pre-training of the model

The training of the neural network is conducted by minimizing the distance
between the error distribution (given by the error model) and the variational
distribution $q_{\theta}$. Here we adopt the Kullback-Leibler divergence.

$\hat{\theta}=\mathop{\arg\min}_{\theta}D_{\mathbf{KL}}\left[P(\alpha,\beta,\gamma)\
||\ q_{\theta}(\alpha,\beta,\gamma)\right]$ (24)

Here we use the forward KL divergence because we always assume that we have
$N$ samples $\\{\alpha,\beta,\gamma\\}\sim P(\alpha,\beta,\gamma)$ of the
errors, which can be obtained by sampling the error model or collected from
experiments. Then the loss function can be evaluated as

$\displaystyle\hat{\theta}$
$\displaystyle=\mathop{\arg\min}_{\theta}D_{\mathbf{KL}}\left[P(\alpha,\beta,\gamma)\
||\ q_{\theta}(\alpha,\beta,\gamma)\right]$
$\displaystyle=\mathop{\arg\min}_{\theta}\sum_{\alpha,\beta,\gamma}P(\alpha,\beta,\gamma)\left[\log
P(\alpha,\beta,\gamma)-\log q_{\theta}(\alpha,\beta,\gamma)\right]$
$\displaystyle=\mathop{\arg\min}_{\theta}\frac{1}{N}\sum_{\\{\alpha,\beta,\gamma\\}\sim
P(\alpha,\beta,\gamma)}\left[\log P(\alpha,\beta,\gamma)-\log
q_{\theta}(\alpha,\beta,\gamma)\right]$
$\displaystyle=\mathop{\arg\min}_{\theta}-\sum_{\\{\alpha,\beta,\gamma\\}\sim
P(\alpha,\beta,\gamma)}\log q_{\theta}(\alpha,\beta,\gamma)$
$\displaystyle=\mathop{\arg\min}_{\theta}F_{\mathrm{NLL}},$ (25)

where $F_{\mathrm{NLL}}=-\sum_{\\{\alpha,\beta,\gamma\\}\sim
P(\alpha,\beta,\gamma)}\log q_{\theta}(\alpha,\beta,\gamma)$ is the so-called
negative log-likelihood loss function or the cross-entropy loss.

In this work, we consider error models that are easy to sample. For example,
for the depolarizing model, the errors on each qubit are generated
independently. We also consider the correlated noise where the errors are
generated according to some pairwise correlations. In this case, we can adopt
the Metropolis-Hasting algorithm to sample the error. Once an error operator
is sampled, according to the isomorphic mappings between the Pauli group and
its ELS representation (i.e. with Eq. 15, Eq. 16, and Eq. 17), a
$\\{\alpha,\beta,\gamma\\}$ configuration is generated as training data for
training the model for all syndromes. So we call it Pre-training, it learns a
joint distribution for all syndromes, and also keeps the conditional
probabilities $q(\beta|\gamma)$ for all syndromes.

So the decoding is very fast because one can evaluate the conditional
probabilities for each $\beta_{i}$ variable by the forward passes of the
neural network especially using GPUs. Here we describe it in detail using the
example of Fig. 3. First, a syndrome $\gamma_{1},\gamma_{2}$ is sent as an
input to the Transformer, after one forward pass of the Transformer
$F_{\beta_{1}}(\gamma_{2},\gamma_{1})$, we compute a conditional probability
$q(\beta_{1}|\gamma_{2},\gamma_{1})$ using e.g. Eq. 21 and sample a
configuration of $\beta_{1}$ which maximize the conditional probability with

$\hat{\beta}_{1}=\mathop{\arg\max}_{\beta_{1}}q(\beta_{1}|\gamma_{2},\gamma_{1}).$
(26)

Then, we send $\beta_{1}$ and $\gamma$ as input to the transformer, compute a
conditional probability $q(\beta_{2}|\beta_{1},\gamma_{2},\gamma_{1})$, then
sample a $\beta_{2}$ configuration according to the conditional probability

$\hat{\beta}_{2}=\mathop{\arg\max}_{\beta_{2}}q(\beta_{2}|\beta_{1},\gamma_{2},\gamma_{1}).$
(27)

We call the method of sample $\beta$ variables one by one the generative MLD
decoding. In Fig. 4, we compare the logical error rate given by the exact MLD
decoding which enumerates all possible $4^{k}$ logical operators, and the
generative MLD and we can see that on the surface code and we can see that the
generative results are almost identical to the exact MLD results while the
computational complexity has decreased from $4^{k}$ to $2k$ for $2k$
conditional probabilities.

Figure 4: Comparison between the exact maximum likelihood decoding (MLD) and
an exact MLD (solid lines) and generative MLD (stars). The coset probabilities
$p(\mathbf{\beta}|\gamma)$ and the conditional coset probabilities
$p(\beta_{i}|\beta_{j<i},\gamma)$ are calculated through contraction of tensor
networks. We do the lattice surgeries on d=5 Surface Code. The numbers of
logical qubits are [1, 3, 5]. Each data point is averaged over $10000$ random
syndrome instances.

We remark that the decoding for multiple syndromes can be done simultaneously
because the Transformer can accept a batch of syndromes as an input and
process the batched forward pass efficiently, especially using GPUs.

##  Appendix G Refinement

The pre-training learns a joint distribution $q_{\theta}(\alpha,\beta,\gamma)$
for all $2^{m}$ $\gamma$ configurations and offers a fast forward pass for
decoding. For a given syndrome, we can spend more computational cost to
further enhance the accuracy of decoding, we term it as refinement. A
straightforward way for the refinement of the transformer is minimizing the
distance between the conditional distribution and the true conditional
distribution. However, since it is not possible to obtain samples of the
conditional distribution $p(\alpha,\beta|\gamma)$ given a syndrome directly
from the error model due to the lack of the normalization factor, we can not
directly refine the variational distribution given the syndrome using the
forward KL divergence
$D_{\mathrm{K}L}(p(\alpha,\beta|\gamma)\|q(\alpha,\beta|\gamma)$. Instead, we
can do the refinement by minimizing the reverse KL

$D_{\mathrm{K}L}(q(\alpha,\beta|\gamma)\|p(\alpha,\beta|\gamma)),$

because the variational distribution $q(\alpha,\beta|\gamma)$ is always sample
able. However, using the reverse KL divergence requires computing the
gradients using the samples of the variational distribution with the REINFORCE
method [28] (also known as policy gradients), analogous to reinforcement
learning. The procedure is quite similar to the variational autoregressive
neural networks for minimizing the variational free energy for statistical
mechanics problems [29].

However, minimizing the backward KL is computationally expansive. In this
work, we propose another way for the refinement of qecGPT, which we call the
generative refinement. The idea is to explicitly compute the summation of the
stabilizer configurations using the $\alpha$ configurations sampled from the
variational distribution (taking advantage of the generative models), given a
syndrome configuration $\gamma$. That is the unbiased version of the joint
probability $p(\beta,\gamma)$ can be computed using samples of stabilizer
configurations that are reweighted.

$\displaystyle p(\beta,\gamma)$
$\displaystyle=\sum_{\alpha}q_{\theta}(\alpha,\beta,\gamma)\frac{p(\alpha,\beta,\gamma)}{q_{\theta}(\alpha,\beta,\gamma)}\approx\frac{1}{N}\sum_{\alpha\sim
q}\frac{p(\alpha,\beta,\gamma)}{q_{\theta}(\alpha,\beta,\gamma)}.$

On the R.H.S. of the last equation, we use the samples of the variational
distribution and the reweighting to compute an unbiased estimate of the joint
distribution. The samples of $\alpha$ configurations can be computed using
conditional probabilities, but different from the way we obtained $\beta$
configurations. Again we use Fig. 3 as a simple example. Suppose we are given
a syndrome configuration $\gamma_{1},\gamma_{2}$ and we want to estimate the
joint probability of a $\beta$ configuration $\beta_{1},\beta_{2}$ and the
syndrome configuration $p(\beta_{2},\beta_{1},\gamma_{2},\gamma_{1})$. We can
send $\beta_{2},\beta_{1},\gamma_{2},\gamma_{1}$ as an input to the
transformer, compute the conditional probability and sample the $\alpha_{1}$
configuration according to this probability distribution

$\alpha_{1}\sim q(\alpha_{1}|\beta_{2},\beta_{1},\gamma_{2},\gamma_{1}).$

Notice that this is different from the sampling procedure to determine the
$\beta$ configuration in the decoding of the pre-trained model, as here we are
sampling from the distribution while in decoding of the pre-trained model
maximizing the conditional probabilities as shown in Eq. (26). After
determining the value of $\alpha_{1}$, we send
$\alpha_{1},\beta_{2},\beta_{1},\gamma_{2},\gamma_{1}$ as an input to the
transformer, compute the conditional probability and sample the $\alpha_{1}$
configuration according to this it

$\alpha_{2}\sim
q(\alpha_{2}|\alpha_{1},\beta_{2},\beta_{1},\gamma_{2},\gamma_{1}).$

The advantage of the reweighting formula is the unbiased estimates, and it
greatly improves the decoding accuracy from the pre-training, but the drawback
is that we have to evaluate $4^{k}$ configurations of $\beta$ variables for
the maximum-likelihood decoding.

##  Appendix H Fast decoding with mismatched training parameters to the error
model

In maximum likelihood decoding, computation of the likelihood of logical
operators requires the parameters of the error model. In the qecGPT, the
parameters of the error model explicitly appear in the (slow) training process
but do not appear in the (fast) decoding process. It would be very efficient
if a transformer trained using a set of parameters of the noise model but used
for decoding under another noise model parameter. In the area of statistical
inference, this is known as inference with mismatched parameters and it is
well known that although not optimal, sometimes the mismatched parameters
already provide very accurate inference results. For example, in the community
detection problem (which can be studied using the statistical inference of the
stochastic block model), It has been shown that the inference using the
parameters at the phase transition always gives good results, and is even
optimal in the sense of the range of detection [22]. Inspired by the results
of [22], we can always train the qecGPT using the parameters of the noise
model at the theoretical phase transition point. To give a concrete example,
for the surface code, the phase transition happens with $n\to\infty$ at a
physical error rate $p\approx 0.189$ [30]. We can train the qecGPT with
$p=0.189$ and decode for codes with other physical error rate values.

In Fig. 5, we have tasted the performance of mismatched decoding of the
surface code with various code distances. The maximum likelihood decoding is
performed using exact tensor network contractions with noise parameter
$p^{\prime}$. We have tested two kinds of $p^{\prime}$ values. The first one
is the matched parameters, where $p^{\prime}$ is set to the true physical
error rate $p$ in the depolarizing noise model; the second one is fixed to
$p^{\prime}=0.189$ which is the threshold (phase transition) of the surface
code with $d=\infty$. From the figure, we can see that the results of the MLD
decoding with the threshold parameter is indistinguishable from the exact MLD
decoding with true parameters.

Figure 5: The performance of exact maximum likelihood decoding with mismatched
parameters on the $d=[3,7,11]$ surface codes under the depolarizing error
model. The decoding is implemented using exact tensor network contractions
with noise parameter $p^{\prime}$. The blue line shows the decoding results
with $p^{\prime}=p$, i.e. with matched error model parameters. The red symbols
show the decoding results with $p^{\prime}=0.189$, which is fixed to the
physical error rate at the threshold (phase transition). Each data point is
averaged over $10000$ random syndrome instances.

###  H.1 Compared with other neural network decoders

Recently, several decoding algorithms based on neural networks have been
proposed. These include algorithms using Boltzmann machines [7], multilinear
neural networks [8, 9, 10, 11], the Long Short-Term Memory (LSTM) neural
networks [12], and the convolutional neural network (CNN) [13, 14]. We notice
that all of the existing neural network decoders are based on supervised
learning. This means that the training of the models requires a training
dataset with labels prepared using another decoder. The labels are either the
operators that are computed using the minimum weight perfect matching
algorithm, or the correct type of logical operators that are computed using a
maximum likelihood decoding algorithm. For detailed introductions to the
neural network decoders we refer to [31].

In this sense, our generative decoders are much different from the existing
neural network decoders. Among many differences, the most crucial difference
is that our approach uses unsupervised learning rather than supervised
learning. It models the joint distribution of errors using neural networks,
rather than learning the probability of labels. In other words, for training
the neural networks, we do not need to prepare any labeled data. Instead, we
directly draw unlabelled samples from the noise model as training data.
Moreover, the existing neural decoders belong to discriminative learning while
our approach belongs to generative learning which generates the maximum likely
logical operators variable by variable, in analogous to generating a sentence
word by word. Another significant difference between our approach and the
existing neural network decoder is that our generative decoder can generate
maximum likely logical operators variable by variable, hence is capable of
decoding with a large number of logical qubits. On the opposite, the maximum
likelihood decoders based on supervised learning require labeled data with
$4^{d}$ different kinds of labels, which is intractable. For example, with
$k=7$, our generative decoder can generate logical operators with
computational complexity $O(14)$, while the neural network decoders based on
the generative learning require $4^{7}=16384$ different types of labels and
hence is intractable.

##  Appendix I Additional numerical results

In this section, we provide numerical results in addition to the results we
have shown in the main text.

####  Rotated Surface Codes—

Here we benchmark qecGPT on the Rotated Surface Codes [32]. This type of code
has the same threshold as the surface Code but with only $d^{2}$ physical
qubits for encoding a single logical qubit, rather than $n=d^{2}+(d-1)^{2}$
qubits in the surface code. Recently this code is frequently used in quantum
hardware experiments [3]. We compared the performance of qecGPT with MWPM in
Fig. 6 on the rotated surface code with different code distances. It can be
seen that qecGPT always provides more accurate results (with lower logical
error rates) than MWPM.

Figure 6: Comparison of decoding performance of various algorithms on the
rotated surface code with $d=3$(a), $d=5$(b), and $d=7$(c). Each data point is
averaged over $10,000$ random instances. MWPM means minimum weight perfect
matching algorithm, and the Exact-MLD is implemented by summing all possible
stabilizer configurations using exact tensor network contractions.

####  With $k>1$ logical qubits—

An advantage of qecGPT is that it can decode with $k\gg 1$ logical qubits
efficiently via generative capability. In maximum likelihood decoding, one
usually needs to compare the likelihood of all $4^{k}$ logical operators with
$k$ logical qubits, which is intractable for a large $k$. Instead, the
complexity of qecGPT is only $O(2k)$ because it generates logical
configurations variable by variable. To evaluate the performance of qecGPT
with $k>1$ logical qubits, we did lattice surgery on $d=3$ surface code to
increase the number of logical qubits $k$. In detail, we have removed 2 (4, 6)
stabilizers randomly and increased the number of logical qubits to 3 (5, 7)
respectively. This increases the number of logical operators to $64$ with
$k=2$, $1024$ with $k=5$, and $16384$ with $k=7$. During the training of the
transformer, we always fix the physical error rate to $p^{\prime}=0.15$ and
use the model to decode with other physical error rates $p$. From Fig. 7 we
can see that qecGPT outperforms minimum weight perfect matching on both
surface code and toric code [4] with multiple logical qubits.

Figure 7: Comparison of decoding performance with $k>1$ logical qubits.(left)
distance-3 Surface Codes with random lattice surgeries. The number of physical
qubits is 26 and the numbers of logical qubits are k=3 (black), k=5 (blue),
and k=7 (red) respectively. (right) Toric Code with code distance $d=3$ and
$k=2$ logical qubits. In the figures, MWPM is the minimum weight perfect
matching [5], and BPOSD is belief propagation augmented by ordered statistics
decoding (BPOSD) [23] algorithm. Exact MLD is the maximum likelihood decoding
which sums all possible stabilizer configurations using exact tensor network
contractions.

####  Correlated noise—

Although simple error models e.g. depolarizing models assume that errors are
independent on each qubit, in practical quantum hardware there inevitably
exist correlations between errors on different qubits. In our approach, we
only need samples from the error models for training qecGPT, so we can decode
error-correcting codes with correlated noise in exactly the same way as we
described and evaluated with the independent noise models. In this section, we
propose a simple noise model to evaluate our approach for the correlated
noise. For a $[n,k,d]$ quantum code, we generate an Ising model on $degree=4$
regular random graph with $n$ spins. The couplings $J_{ij}$ of the Ising model
are sampled from a uniform distribution $\mathcal{U}(0,1)$. And for breaking
the $Z_{2}$ symmetry, a small external field $h=0.3$ is added to Hamiltonian
as Eq. 28.

$H=-\beta\sum_{<ij>}J_{ij}s_{i}s_{j}-h\sum_{i}s_{i}$ (28)

Then we draw samples from the Boltzmann distribution using the Metropolis-
Hasting algorithm. The error that occurred is determined by a sample $s$. If
$s_{i}=1$, an identity $I$ acts on $i$th physical qubit. Else, There is a
Pauli error $\\{X,Y,Z\\}$ with the same probabilities. We can see that the
error configurations are mapped from the configuration of Ising models which
are generated from the Boltzmann distribution with long-range correlations.
With $\beta=0$, all the configurations are drawn randomly, each spin has a
probability of $0.5$ to take $1$ or $-1$, which is analogous to the
depolarizing error model with a high physical error rate. With $\beta>0$ the
correlations between the error configurations are long-range and non-trivial.
With $\beta=\infty$, the ground state of the Ising model is a ferromagnetic
configuration, so there are almost no Pauli error appears, corresponding to a
low physical error rate. In Fig. 8 we compare the performance of qecGPT with
the minimum weight perfect matching algorithm which determines the weights
using the marginal probabilities calculated using the samples. In the figure,
with $\beta=0$, the error model is purely random so the logical error rate of
both MWPM and qecGPT is high. With $\beta$ large, the physical error rate is
small and the decoding is very easy we also see that both MWPM and qecGPT give
a very low logical error rate. In the middle when the physical error rate is
moderate, we can see that qecGPT significantly outperforms the minimum weight
perfect matching algorithm.

Figure 8: Decoding of the surface code with $d=5$ with correlated errors
described in the text. The weights of MWPM are determined using the marginal
probabilities calculated using samples. Each data point is averaged over
$10000$ instances.

###  I.1 Parameters of the neural networks

In the numerical experiments, all neural networks were trained on a single
NVIDIA A100 GPU. To ensure that the distribution of the Transformers has the
autoregressive property, an upper-triangle mask is added to the attention
block as illustrated in Fig. 3. We note that in the Softmax layer of the
Transformers, we have added a mask matrix with diagonal elements set to $1$.

$\left(\begin{matrix}1&$-$\infty&\cdots&$-$\infty\\\
\vdots&\ddots&\ddots&\vdots\\\ \vdots&&\ddots&$-$\infty\\\
1&\cdots&\cdots&1\\\ \end{matrix}\right)$ (29)

In this setting, the autoregressive property is not satisfied because the
$i$th output depends on the $i$th input. The solution is to introduce a
virtual variable $x_{0}$ and modify the order of input variables as

$\mathbf{X}(x_{1},\cdots,x_{n})\rightarrow\mathbf{X}(x_{0},\cdots,x_{n-1}),$
(30)

and $x_{0}$ always equals 1. The hyperparameters for training the Transformers
are detailed in Tab. 1. Where $D$ is the dimension of the model, $N_{h}$ is
the number of heads of multi-head attention, $N_{l}$ is the number of encoder
layers, and $D_{f}$ is the dimension of feed-forward layers.

| BATCH | EPOCH | LR | $D$ | $N_{h}$ | $N_{l}$ | $D_{f}$
---|---|---|---|---|---|---|---
Sur3 | $10^{4}$ | $10^{5}$ | $10^{-3}$ | 256 | 4 | 2 | 256
Sur5 | $10^{4}$ | $2\times 10^{5}$ | $10^{-3}$ | 256 | 4 | 3 | 256
RSur3 | $10^{4}$ | $10^{5}$ | $10^{-3}$ | 256 | 4 | 2 | 256
RSur5 | $10^{4}$ | $2\times 10^{5}$ | $10^{-3}$ | 256 | 4 | 3 | 256
RSur7 | $10^{4}$ | $3\times 10^{5}$ | $10^{-3}$ | 512 | 4 | 2 | 512
3DSur2 | $10^{4}$ | $10^{5}$ | $10^{-3}$ | 256 | 4 | 2 | 256
Tor3 | $10^{4}$ | $10^{5}$ | $10^{-3}$ | 256 | 4 | 2 | 256
N13k3 | $10^{4}$ | $10^{5}$ | $10^{-3}$ | 256 | 4 | 2 | 256
N13k5 | $10^{4}$ | $10^{5}$ | $10^{-3}$ | 256 | 4 | 2 | 256
N13k7 | $10^{4}$ | $10^{5}$ | $10^{-3}$ | 256 | 4 | 2 | 256
Table 1: Parameters of qecGPT and hyperparameters in the training.