<figure><img src="image_1.png" /></figure>

HAL
open science

# Optimization of Hypergraph Product Codes (HGP) for the quantum erasure channel via Reinforcement Learning

Bruno Costa Alves Freire

► To cite this version:

Bruno Costa Alves Freire. Optimization of Hypergraph Product Codes (HGP) for the quantum erasure channel via Reinforcement Learning. Machine Learning [stat.ML]. 2024. hal-04829463v2

HAL Id: hal-04829463

https://inria.hal.science/hal-04829463v2

Submitted on 12 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH INTERNSHIP REPORT

Optimization of Hypergraph Product Codes
(HGP) for the quantum erasure channel via
Reinforcement Learning

2 Mai 2024 - 30 Septembre 2024

Bruno Costa Alves Freire

École Normale Supérieure Paris-Saclay, Master MVA

Tuteur de Stage : Prof. Dr. Anthony LEVERRIER

Enseignant Référent : Prof. Dr. Pierre GAILLARD

Laboratoire : Équipe COSMIQ, Inria Paris
Abstract

In the field of Quantum Computing, it is widely understood that the full potential of quantum
computers can only be achieved through the implementation of Quantum Error Correction
procedures. A particularly interesting class of quantum codes is that of the quantum Low-
Density Parity-Check (qLDPC) codes, as their properties allow for the implementation of fault-
tolerant quantum computers with constant overhead. However, finding good qLDPC codes
remains a challenge. Such codes can be produced through a particular construction known
as HyperGraph Product (HGP) codes, which are constructed from two classical LDPC codes,
yielding quantum codes with constant rate and relatively good robustness to noise. By adopting
HGP codes, the problem of optimizing a quantum code is reduced to that of optimizing two
classical codes. Here, we formulate the problem of optimizing quantum codes as a search
problem in a space of classical Tanner graphs, to be traversed by means of local deformations to
their connectivity structure. By applying search-based algorithms such as Simulated Annealing
and a Reinforcement Learning algorithm known as Projective Simulation, we optimize quantum
codes for the quantum erasure channel, obtaining performances up to one order of magnitude
better than classical methods like the Progressive-Edge Growth algorithm.

Résumé

Dans le domaine de l'informatique quantique, il est largement admis que le plein potentiel des
ordinateurs quantiques ne peut être atteint qu'à travers la mise en œuvre de procédures de
correction d'erreurs quantiques. Une classe particulièrement intéressante de codes quantiques
est celle des codes quantiques à faible densité de contrôle de parité (qLDPC), car leurs proprié-
tés permettent la mise en œuvre d'ordinateurs quantiques tolérants aux fautes avec un surcoût
constant. Cependant, trouver de bons codes qLDPC reste un défi. De tels codes peuvent être
produits grâce à une construction particulière connue sous le nom de codes produits d'hyper-
graphes (HGP), construits à partir de deux codes LDPC classiques, produisant ainsi des codes
quantiques à taux constant et relativement robustes face au bruit. En adoptant les codes HGP,
le problème d'optimisation d'un code quantique se réduit à celui de l'optimisation de deux
codes classiques. Ici, nous formulons le problème de l'optimisation des codes quantiques comme
un problème de recherche dans un espace de graphes de Tanner classiques, à traversé par des
déformations locales de leur structure de connectivité. En appliquant des algorithmes basés sur
la recherche tels que le Recuit Simulé et un algorithme d'apprentissage par renforcement connu
sous le nom de Simulation Projective, nous optimisons des codes quantiques pour le canal d'ef-
facement quantique, obtenant des performances jusqu'à un ordre de grandeur supérieur aux
méthodes classiques comme l'algorithme de croissance progressive des arêtes (Progressive-Edge
Growth).
Acknowledgement

The greatest acknowledgements of this work are due to my supervisor, Prof. Dr. Anthony LEVERRIER, for accepting the challenge of dealing with an unfamiliar field of research, that of machine learning methods, all the while providing careful guidance and extremely helpful advice on the conduction of the work.

I thank Nicolas DELFOSSE for accepting to collaborate in this work by sharing his experi-
ence with related works, and by pointing out solutions for technical implementation questions.

I'm also grateful towards Prof. Dr. Jean-Pierre Tillich, for having introduced me to the
world of Information Theory and Error Correcting Codes, as well as serving as a person of
contact between me and my supervisor.

Last but not least, I am grateful to the CLEPS infrastructure from the Inria of Paris for
providing resources and support.
# Contents

<table><tr><td>Acknowledgement</td><td>1</td></tr><tr><td>Introduction</td><td>3</td></tr><tr><td>1 Quantum Error Correction</td><td>6</td></tr><tr><td>1.1 Elements from Classical Error Correction</td><td>6</td></tr><tr><td>1.1.1 Linear codes</td><td>7</td></tr><tr><td>1.2 Basics of Quantum Information</td><td>9</td></tr><tr><td>1.2.1 Qubits</td><td>9</td></tr><tr><td>1.2.2 Quantum measurements</td><td>11</td></tr><tr><td>1.3 Core concepts of Quantum Error Correction</td><td>12</td></tr><tr><td>1.3.1 Stabilizer formalism</td><td>12</td></tr><tr><td>1.3.2 CSS codes</td><td>15</td></tr><tr><td>1.3.3 The quantum erasure channel</td><td>16</td></tr><tr><td>1.4 qLDPC and HGP codes</td><td>18</td></tr><tr><td>2 Optimizing QEC codes</td><td>21</td></tr><tr><td>2.1 Modeling the problem in RL framework</td><td>21</td></tr><tr><td>2.1.1 State space</td><td>22</td></tr><tr><td>2.1.2 Action space</td><td>23</td></tr><tr><td>2.1.3 Cost function</td><td>24</td></tr><tr><td>2.2 Projective Simulation</td><td>26</td></tr><tr><td>2.2.1 Implementation details</td><td>27</td></tr><tr><td>2.3 Simulated Annealing and plain exploration</td><td>29</td></tr><tr><td>3 Results and Discussion</td><td>30</td></tr><tr><td>3.1 Plain Exploration</td><td>30</td></tr><tr><td>3.2 Simulated Annealing</td><td>32</td></tr><tr><td>3.2.1 Plain acceptance probability function</td><td>32</td></tr><tr><td>3.2.2 Log-scale acceptance probability function</td><td>33</td></tr><tr><td>3.3 Projective Simulation</td><td>34</td></tr><tr><td>3.4 Critique of the methods, implementations and results</td><td>36</td></tr><tr><td>Conclusion</td><td>40</td></tr><tr><td>References</td><td>41</td></tr></table>
# Introduction

In recent years, the field of Quantum Computing has received increasing attention from both industry and academia. Many physical platforms for implementing quantum computers have been proposed and are under active research to improve the underlying technologies and to probe their inherent limitations. While some researchers and companies investigate the capabilities and potential uses for devices in what is called the *Noisy Intermediate-Scale Quantum (NISQ)* era, it is widely believed that the majority of the applications of interest for the quantum computing paradigm can only be fully exploited through *Fault-Tolerant Quantum Computing (FTQC)*.

The distinction between the two lies primarily in the level of reliability of the quantum devices, measured essentially by the error rate of the system, roughly speaking, how long the system can operate, or how many quantum operations it can undergo without being compromised by some sort of noise-induced error. In practical terms, the existing approaches that try to achieve FTQC rely either on inherent physical protections against noise, such as topological quantum computers [1], bosonic qubits [2, 3], to name a few, or, more broadly, on Quantum Error Correction (QEC), which roughly consists on encoding a number of logical qubits into a larger number of physical qubits, in a way that trades redundancy for robustness against noise. One could make an analogy with hardware-based versus software-based solutions.

Currently there are many different methods of implementing QEC, and there is active research both on the formulation of new families of QEC codes and on the design of decoding algorithms. On the side of the decoders, there is a wide literature on algorithms focusing on approximate, but fast methods [4], for different quantum channels and for different code families [5, 6, 7]. There is also plenty of research applying machine learning techniques to the development of decoders, such as neural networks and reinforcement learning based decoders [8, 9, 10]. However, in the end of the day, the error suppression capabilities of a QEC scheme are still limited by the encoder side. In this sense, this work seeks to optimize the encoding part of this endeavor.

A particularly interesting class of QEC codes is that of the quantum Low-Density Parity-Check (qLDPC) codes [11]. Like their classical analogues, these codes are characterized by parity-check matrices with low densities, meaning that the measurement of each syndrome component acts only on a small fraction of the qubits of the system, and each qubit is involved in a small fraction of the syndrome component measurements. With each qubit being involved in a constant number of parity-check equations, this property enables a syndrome measurement circuit of constant depth, and the efficiency in this step together with a fast classical decoder allow for FTQC with constant overhead [12]. In other words, it enables the suppression of errors up to a desired level with only a linear increase in the number of qubits. However, in the quantum case, qLDPC codes are still not very well understood, in particular because it is difficult to assess their performance in realistic noise models, due to the lack of efficient optimal decoders.

The current approach to optimizing qLDPC codes consists on considering some code construction such as the Hypergraph Product (HGP) codes [13]. These codes are interesting, for instance, due to their compatibility with some recent quantum hardware architectures based
on neutral atoms [14]. Moreover, HGP codes also profit from good asymptotic properties, namely an asymptotically constant rate of encoding, and a code distance proportional to the square root of the number of qubits, this last property being a measure of robustness against adversarial noise. Recently, code constructions yielding even better properties have been discovered, such as lifted product codes [15] and quantum Tanner codes [16], which have both linear dimension (constant rate) and linear distance, and generalize the HGP construction. For simplicity reasons, this work explores only HGP codes.

HGP codes are produced from two classical codes by taking the cartesian product of their Tanner graphs. When the input classical codes are LDPC, the result is a qLDPC code. The problem is then reduced to optimizing the two classical factor codes. For this, the current strategies are mostly ad-hoc, such as trying to reduce the number of short cycles, or removing small trapping sets, which are both elements that negatively impact the quality of decoding algorithms in the classical setting. An example of a standard technique in this sense is the Progressive Edge-Growth algorithm [17], which aims at producing classical LDPC codes whose Tanner graphs have a large girth, which is the length of the shortest cycle in the graph.

## Review of related applied ML works

In such a context, it is interesting to investigate the applicability of machine learning techniques to the discovery and optimization of qLDPC codes. In [18], the authors apply Proximal Policy Optimization (PPO) to the task of discovering quantum codes and their encoding circuits. The state space is modeled by Pauli strings describing generators of the stabilizer group (see section 1.3.1 for the definition), and the actions consist on adding a next gate to the encoding circuit. The reward in their case is based on a verification of the Knill-Laflamme conditions, which is roughly whether or not a given error pattern can be detected by the error correction code. This is the source of their main limitation, since the verification of these conditions scales exponentially with the number of qubits.

Other works have also explored RL techniques to the task of optimizing and discovering QEC codes. In [19], for instance, the authors use a different framework of QEC codes, known as Quantum LEGOs, and apply Masked Proximal Policy Optimization (MPPO) as their RL algorithm. Their rewards are based on adversarial distance (the minimum distance of the code) and on the logical error rate under a biased noise model. Because they compute the logical error rate exactly, by means of the quantum weight enumerator polynomials, their approach also scales exponentially, limiting their method to relatively small numbers of qubits. [20], on the other hand, works with Tensor Network codes, and utilizes an RL framework known as Projective Simulation (PS), which is a physics-inspired AI framework introduced in [21] and [22]. Although less known in the RL community, the PS framework is simple and flexible. In [20], the state is again given by the set of generators of the stabilizer group describing the code, but the actions here consist on tensor-network operations, such as node contractions. The reward is based on the minimum distance of the codes, which are pre-computed for codes up to some small size.

In our work, an approach similar to that of [23] is adopted. In that paper, the authors also employ Projective Simulation, in a setting where the state space is formed by a family
of surface codes (see [24] for definition), and the action space consists of local deformations of such codes. The reward model is based on the evaluation of the code's performance under a particular quantum channel, which is enabled by a linear-time benchmarking tool developed in [25], which is specific to surface codes and relies on a simplification of the noise model. By considering a simplified noise model, in this case the quantum erasure channel, one can simulate the performance of a qLDPC code with reasonable efficiency, thus allowing for the use of a similar approach.

In view of the related works, the adopted approach for this internship was to employ Pro-
jective Simulation, tailoring the state and action spaces to the context of HGP codes. The
reward model was based one the Monte Carlo estimation of the logical error rate of the code
under the quantum erasure channel. This choice can be justified by analogy with the classical
case, where the study of codes for the erasure channel served as a proxy of their performance
in other, more realistic noise models. In addition to exploring PS, we develop a comparison of
this framework with other optimization techniques such as Simulated Annealing.

The report is structured as follows: in section 1, some background on QEC is provided. In section 2, the optimization methods explored are reviewed and the QEC optimization problem is tailored to the adopted framework. Section 3 presents and discusses the results and points out further directions of investigation. All the code developed to implement the modeling of the optimization problem, the optimization algorithm and the running of experiments and generation of plots is gathered in the project's github repository.
# 1 Quantum Error Correction

The aim of this section is to provide some background on Quantum Error Correction for readers unfamiliar with the topic. Fine-grained details are not essential for the understanding of the project, but I opted to present the subject in a rather self-contained, motivated way. Most of the material exposed here is contained in [26]. For a more in-depth treatment of the topics in this section, consult [27].

## 1.1 Elements from Classical Error Correction

Let us introduce some important concepts from classical coding theory that will be mentioned quite often throughout this report, and which will have their quantum analogues showing up as well. In channel coding theory, we have two communicating entities, the sender and the receiver, who must communicate through a channel, which is affected by noise, and can thus corrupt the messages. The sender wants to protect a message consisting of $k$ bits, $u \in \mathbb{F}_2^k$ from the channel's noise, in a way that the receiver can recover the original message with high probability, even arbitrarily close to 1. To do so, both sender and receiver pick a codebook $C \subset \mathbb{F}_2^n$, consisting on $n$-bit codewords, and the sender maps the message $u$ into some codeword $x \in C$ through some encoding function $\phi: \mathbb{F}_2^k \to C$. Since the codebook is only a subset of all the possible messages, the receiver can detect the occurrence of errors whenever the received message $y \in \mathbb{F}_2^n$ is not an element of $C$. The amount of redundancy in the coding scheme is measured by the rate $k/n$.

Speaking of errors, there are two particularly important noise models in classical coding theory, namely the *binary symmetric channel* (BSC), and the *erasure channel*. The BSC can be modeled as a stochastic mapping $x \mapsto y = x + e$, where $e \sim \mathcal{B}(p)^{\otimes n}$ is a random vector in $\mathbb{F}_2^n$ whose entries are i.i.d. bits following a Bernoulli distribution. I.e., the BSC flips each bit of the message independently with probability $p$. The erasure channel, on the other hand, doesn't flip bits but deranges them in such a way that the receiver cannot tell whether it is a 1 or a 0. The modeling of an erasure channel consists on the sampling of an erasure vector $\epsilon \sim \mathcal{B}(p)^{\otimes n}$, and the received message is given by $y \in (\mathbb{F}_2 \cup \{?\})^n$, where $y_i = ? \Leftrightarrow \epsilon_i = 1$. In other words, the bits in the *support* of the erasure vector are converted into wildcards, while the remaining bits are left untouched. See Figure 1 for schematics.

<figure><img src="image_3.png" /><figcaption>(a) Binary symmetric channel with flip rate $p$.</figcaption></figure>

<caption>(b) Erasure channel with erasure rate $p$.</caption>

**Figure 1:** Representation of the most common classical noise models. The diagrams represent the action of the channel on each bit independently.
# 1. QUANTUM ERROR CORRECTION

While both channels can be considered realistic models depending on the communication scenario, it is worth pointing out that the BSC is usually the target model when designing error correction codes, with the erasure channel usually serving as a simplified model that helps analyzing algorithms and understanding the properties of the codes. When LDPC codes were first introduced in [28], they were first analyzed and optimized for the erasure channel, and it served as a proxy to the codes' performance in more realistic channels.

Focusing on error correction in the BSC, the most basic approach to correct such errors is to consider the *closest* codeword to $y$ in the Hamming metric,

$$ \hat{x} \in \arg \min \{|y - w| : w \in \mathcal{C}\}. $$

Of course, such approach is successful provided that this *closest* codeword is unique, and that the error weight $|e|$ is sufficiently small. After all, there is always a nonzero probability that the channel perfectly corrupts one codeword into another, frustrating the chances of the error correction procedure. Knowing that errors with high weight are less likely for a reasonably small $p$, one aims to design codes where the codewords are far apart from each other. This leads to the definition of another important parameter of a code, which is the *minimum distance*

$$ d := \min\{|x - x'| : x \neq x' \in \mathcal{C}\}, $$

which is a measure of robustness of the code against adversarial noise. A code encoding $k$ logical bits into $n$ physical bits, with distance $d$, is called a $[n, k, d]$ code. When the distance is unknown, we denote it as just $[n, k]$.

## 1.1.1 Linear codes

In general one wants to find codes with the highest possible rate and distance. While Shannon's channel coding theorem [29] points out to the existence of random codes achieving the best trade-off between redundancy and robustness, random codes are inconvenient for algorithmic purposes, particularly the decoding part. A class of codes of great interest is that of the *linear codes*. A code is linear when its codebook $\mathcal{C}$ is a linear subspace of $\mathbb{F}_2^n$. In this context, $k$ is called the dimension of the code, and the number of logical bits (as opposed to the number of physical bits, $n$). The distance $d$ is then given by the smallest weight of a nonzero codeword.

Linear codes are convenient for both encoding and decoding steps. In order to encode a message $u \in \mathbb{F}_2^k$, one considers a *generator matrix* $G \in \mathbb{F}_2^{n \times k}$, satisfying $\mathrm{im} G = \mathcal{C}$, i.e. its columns form a basis of the code, and the encoded message is given simply by $x = Gu$. For the decoding step, one considers a dual representation of the code $\mathcal{C}$, by means of a *parity-check matrix* $H \in \mathbb{F}_2^{m \times n}$, satisfying $\mathrm{ker} H = \mathcal{C}$. With this parity-check matrix, the receiver can compute the *syndrome* of the message $y$, $\sigma(y) = Hy \in \mathbb{F}_2^m$, which tells whether the message belongs to the code or not, i.e., $\sigma(y) = 0 \Leftrightarrow y \in \mathcal{C}$. Given the syndrome $\sigma$, the decoder can solve a linear system $H\hat{e} = \sigma (= H(x + e) = He)$ to infer the error and correct it. The equations forming this system are called *parity-check equations* (because they are in $\mathbb{F}_2$), hence the name of $H$. Notice though that this system is underdetermined, since $\mathrm{ker} H = \mathrm{im} G$, i.e. $\mathrm{rank} H = n - k$, so the error $\hat{e}$ still needs to be inferred somehow.
The syndrome extraction trades the problem of inferring the message $x$ to that of inferring the error $e$, and is especially important in the quantum setting, where one cannot measure the quantum state directly. At this point, one natural decoding strategy consists on the *minimum weight decoder*, which infers the error as

$$ \hat{e} := \arg \min \{|e| \in \mathbb{F}_2^n : \sigma(e) = \sigma\}. $$

Linear codes can be represented either by their generator matrix $G$ or by their parity-check matrix $H$, in the sense that each one of those matrices completely characterizes the linear code. Since both are connected by the relation $\mathrm{im} G = \mathcal{C} = \ker H$, one can find each representation from the other, as stated by the following lemma.

**Lemma 1.1** ($\eta-\gamma$ lemma). $\forall G \in \mathbb{F}_2^{n \times k}, \exists \eta(G) \in \mathbb{F}_2^{m \times n}$ such that $\ker \eta(G) = \mathrm{im} G$ and $m = n - \mathrm{rank} G$. Conversely, $\forall H \in \mathbb{F}_2^{m \times n}, \exists \gamma(H) \in \mathbb{F}_2^{n \times k}$ such that $\mathrm{im} \gamma(H) = \ker H$ and $k = n - \mathrm{rank} H$.

The precise algorithm employed to do this conversion will be detailed in section 2.

Another useful representation of a linear code is through its *Tanner graph* $\mathcal{T}(C)$, which is defined as a bipartite graph with a set of $m$ check nodes $C$, and a set of $n$ bit nodes $V$, whose connectivity structure is given by $H$, i.e., the edge $(c_i, v_j)$ belongs to $\mathcal{T}(C)$ if and only if $H_{ij} = 1$. In other words, the Tanner graph $\mathcal{T}(C)$ tells us which bits are acted on by which checks. The $H$ matrix is thus the *biadjacency matrix* of the bipartite graph $\mathcal{T}(C)$. Figure 2 showcases an example of a linear code in its many possible representations.

$$ H = \begin{pmatrix} 1 & 1 & 1 & 0 & 1 & 0 & 0 \\ 0 & 1 & 1 & 1 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 & 1 & 1 & 1 \end{pmatrix} $$

im $G = \mathcal{C} = \ker H$

$$ G = \begin{pmatrix} 1 & 1 & 1 & 0 \\ 1 & 1 & 0 & 1 \\ 0 & 1 & 1 & 1 \\ 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} $$

<figure><img src="image_4.png" /><figcaption>(b) Tanner graph $\mathcal{T}(C)$ of the linear code $\mathcal{C}$.</figcaption></figure>

Figure 2: The multiple different representations of the [7, 4, 3] Hamming code, via its parity-check matrix, its generator matrix and its Tanner graph.
# 1. QUANTUM ERROR CORRECTION

Another important class of linear codes is that of the low-density parity-check (LDPC) codes, introduced by [28]. In general, LDPC codes are codes for which the parity-check matrix $H$ is sparse, but one can also specify more precise LDPC codes, such as $(l, r)$-LDPC codes, which are characterized by the *left* and *right* degrees of their Tanner graphs. Specifically, when a code or a family of codes has their left degrees (i.e., the degrees of the check nodes) and right degrees (bit nodes) bounded respectively by $l$ and $r$, they are said to be $(l, r)$-LDPC.

## 1.2 Basics of Quantum Information

In quantum mechanics, the state of a system is described as an element of a complex Hilbert space, $\mathcal{H}$, called its *wave function* $|\psi\rangle$. Composite systems are modeled by means of tensor products of the Hilbert spaces describing each component system, which gives rise to the phenomenon of *entanglement*, which is when a the state of a composite system cannot be described by the tensor product of states in its component systems, but rather by linear combinations of product states. Quantum states should be normalized, i.e., $|||\psi\rangle|| = \sqrt{\langle\psi|\psi\rangle} = 1$, and moreover, their *global phase* has no physical meaning, that is, for any $z \in \mathbb{S}^1 := \{z \in \mathbb{C} : |z| = 1\}$, the vectors $|\psi\rangle$ and $z|\psi\rangle$ represent the same physical state.

Another way to represent the physical state of a quantum system that captures this feature, is by its *density operator*, $\rho$, which is given by $\rho = |\psi\rangle\langle\psi|$, for *pure* states. In full generality, a density operator is a self-adjoint ($\rho^\dagger = \rho$), positive semi-definite ($\langle u|\rho u\rangle \ge 0 \forall u \in \mathcal{H}$) bounded linear operator, of trace class ($\sum_k \langle e_k|\rho e_k\rangle < \infty$ for any Hilbert basis $\{e_k\}_{k \subset \mathcal{H}}$), with $\text{Tr}(\rho) = 1$. This representation has the flexibility to describe both *pure* and *mixed* states, the latter being statistical mixtures of pure states. Geometrically, the set of all density operators is a convex set (the convex hull of all pure states), and the pure states are its extreme points.

Observable quantities in quantum systems are defined by hermitian (self-adjoint) bounded operators in the underlying Hilbert space. The most important observable of a system is the Hamiltonian, $H$, which gives the energy of a system and governs its temporal evolution. The dynamics of closed quantum systems follows Schrödinger's equation,

$$i\partial_t |\psi(t)\rangle = H(t) |\psi\rangle \quad \text{or} \quad i\partial_t \rho(t) = [H(t), \rho(t)]$$

where $H(t)$ is the Hamiltonian of the system. This yields a form of temporal evolution called *unitary*, meaning that $|\psi(t)\rangle = U(t)|\psi(0)\rangle$, for $U(t)$ a unitary operator ($U^\dagger U = I = UU^\dagger$). This kind of evolution is reversible and in some sense “deterministic”, when thinking about the wave function or density operator of the system. However, for open quantum systems, the evolution of the state is also affected by *measurements*, which will be discussed in 1.2.2.

### 1.2.1 Qubits

In the context of quantum computing, our physical systems consist of *qubits*, or arrays of qubits. A qubit is just a state of a two dimensional system, i.e., $\mathcal{H} = \mathbb{C}^2 = \text{span}\{|0\rangle, |1\rangle\}$, commonly characterized by a basis of two orthonormal states $|0\rangle$ and $|1\rangle$, in analogy to classical bits. Because of the normalization condition and the global phase invariance, qubits can be represented in a 3D way as points in the *Bloch sphere*, see figure 3.
Arrays of N qubits are states in the space $(\mathbb{C}^2)^{\otimes N} \cong \mathbb{C}^{2N}$, with a standard basis given by $|j\rangle$, $j = 0, \dots, 2^N - 1$, which should be interpreted as a tensor product of states $|0\rangle$ and $|1\rangle$ in each entry of the binary representation of j. E.g., for $N = 4$, $|9\rangle = |1001_2\rangle = |1\rangle \otimes |0\rangle \otimes |0\rangle \otimes |1\rangle = |1\rangle |0\rangle |0\rangle |1\rangle$, omitting $\otimes$ for brevity.

<figure><img src="image_5.png" /><figcaption>Figure 3: Representation of a qubit in the pure state $|\psi\rangle = \cos(\theta/2)|0\rangle + e^{i\phi}\sin(\theta/2)|1\rangle$ in the Bloch sphere. Antipodal states in the sphere represent orthogonal states in $\mathbb{C}^2$.</figcaption></figure>

An important set of operators acting on single qubits is the set of Pauli operators, represented in the computational basis $\{|0\rangle, |1\rangle\}$ as

$$X = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}, \quad Y = \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix}, \quad Z = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}.$$

The pairs of antipodal states $|\pm\rangle := \frac{1}{\sqrt{2}}(|0\rangle \pm |1\rangle)$, $|\pm i\rangle := \frac{1}{\sqrt{2}}(|0\rangle \pm i|1\rangle)$ and $(|0\rangle, |1\rangle)$ in figure 3 are the respective eigenstates of Pauli operators X, Y and Z. The three operators are self-adjoint and also unitary, moreover, every unitary U acting on a single qubit can be written as a linear combination of I, X, Y and Z, satisfying a normalization condition. For multiple qubits, the same result holds for tensor products of Pauli operators. The Pauli operators also profit from a group structure inside the unitary group. The n-qubit Pauli group, that plays a fundamental role in QEC, is defined by

$$\mathcal{P}_n = \{i^r P : r = 0, 1, 2, 3, P \in \{I, X, Y, Z\}^{\otimes n}\}.$$

The group structure of $\mathcal{P}_1$ is summarized by

$$\left\{ \begin{array}{l} X^2 = Y^2 = Z^2 = I \\ XYZ = iI \end{array} \right. .$$
# 1. QUANTUM ERROR CORRECTION

It is worth pointing out that every two elements in the Pauli group either commute or anti-commute, in fact X, Y, Z all anti-commute with each other, and commute with I and with themselves. For QEC, it is also fruitful to consider the phaseless Pauli group, $\tilde{P}_n = P_n/\{\pm 1, \pm i\} \cong K^{\otimes n}$, where K is the Klein group.

As unitary operators acting on qubits, and also as observables with close connections to the computational basis, the Pauli operators play an important role in what concerns the dynamics of qubit systems.

## 1.2.2 Quantum measurements

To complete the description of temporal evolution in qubit systems, we need cover the effect of measurements. To describe measurements mathematically, one considers the framework of projective measurements, which consist on a set of orthogonal projections $\{\Pi_k\}_k$ that act on a state by the mapping $\rho \mapsto \sum_k \Pi_k \rho \Pi_k^\dagger$. This framework can be generalized to positive operator-valued measurements (POVMs) by allowing the operators $\Pi_k$ to be simply positive semi-definite operators. For concreteness, consider a set of orthogonal projectors acting on n qubits, $\{\Pi_k\}_k \subset \mathbb{C}^{2^n \times 2^n}$, satisfying $\sum_k \Pi_k = I$. The process measurement yields a piece of classical information $\mu$, and randomly collapses the quantum state $|\psi\rangle$ into one of the pure states

$$ \left\{ \frac{\Pi_k |\psi\rangle}{\sqrt{\langle\psi|\Pi_k|\psi\rangle}}, \text{ yielding } \mu = k \text{ with probability } \langle\psi|\Pi_k|\psi\rangle \right. . $$

Since this evolution is stochastic, the resulting state is a statistical ensemble of pure states, called a mixed state. A mixed state is basically a convex combination of pure states, and in this case, it is given in terms of the density operators as a sum weighted by the probabilities of each outcome, $|\psi\rangle\langle\psi| = \rho \mapsto \sum_k \Pi_k \rho \Pi_k^\dagger$. As an example, consider the measurement of the observable Z. The associated projection operators are $\Pi_0 = (I+Z)/2$ and $\Pi_1 = (I-Z)/2$, the spectral projectors onto each eigenstate of Z. For a pure state $|\psi\rangle = \alpha|0\rangle + \beta|1\rangle$, the outcome of the measurement is

$$ \left\{ \begin{array}{ll} |\psi'\rangle = |0\rangle, & \text{yielding } \mu = 0 \text{ with probability } |\alpha|^2 \\ |\psi'\rangle = |1\rangle, & \text{yielding } \mu = 1 \text{ with probability } |\beta|^2 \end{array} \right. . $$

More generally, a Pauli measurement associates to each Pauli operator acting on n qubits P, the projections $\Pi_0 = (I+P)/2$ and $\Pi_1 = (I-P)/2$.

To wrap up the dynamics of qubit systems, we introduce the notion of a *quantum channel*, which encompasses all valid evolutions of a system of multiple qubits. It is given by a completely-positive and trace preserving (CPTP) map $\mathcal{E}$, defined by a set of operators $\{E_k\}$ satisfying $\sum_k E_k E_k^\dagger = I$, called *Kraus operators*. The action of $\mathcal{E}$ is given by $\mathcal{E}: \rho \mapsto \sum_k E_k \rho E_k^\dagger$. It's easy to see how this definition entails both measurements and unitary evolutions. A central example of such a channel is the *depolarizing channel*, described by the Kraus map

$$ \mathcal{E}_p : \rho \mapsto (1-p)\rho + \frac{p}{3}X\rho X + \frac{p}{3}Y\rho Y + \frac{p}{3}Z\rho Z. $$
Another important model of quantum channel is the erasure channel

$$ \mathcal{E}_p : \rho \mapsto (1-p)\rho + p|?\rangle\langle?|, $$

which will be discussed in section 1.3.3.

## 1.3 Core concepts of Quantum Error Correction

In classical error correction, one only has to worry about one kind of errors, the bit flips. In the quantum case, however, qubits may suffer from a continuum of unitary errors. Individual qubits may wander around the Bloch sphere, and they may also get entangled with each other. To perform quantum error correction in a “digital” way seems impossible, one might wonder. However, on the foundation of quantum error correction codes lie the aforementioned *Pauli measurements*. The key insight to correct continuous errors in a discrete way is that, by applying Pauli measurements, the state of the system is collapsed, resulting in errors that have the same effect as Pauli operators, rendering the error model discrete again.

### 1.3.1 Stabilizer formalism

A quantum error correction code $\mathcal{Q}$ is defined, in analogy to classical linear codes, as a linear subspace of $\mathcal{H} = \mathbb{C}^{2^N}$ of dimension $2^K$, where $N$ is the number of physical qubits and $K$ is the number of logical qubits. In the stabilizer formalism, introduced by Gottesman in [30], one considers the set of all Pauli operators $S \in \mathcal{P}_N$ that fix each element of the codespace $\mathcal{Q}$:

$$ \mathcal{S} = \{S \in \mathcal{P}_N : S|\psi\rangle = |\psi\rangle, \forall |\psi\rangle \in \mathcal{Q}\} \qquad (1) $$

Notice that the set $\mathcal{S}$ constitutes a subgroup of $\mathcal{P}_N$, as it contains $I$ and is closed under multiplication and inversion (recall that every element of the Pauli group is its own inverse). We call $\mathcal{S}$ the stabilizer group of the code $\mathcal{Q}$. Conversely, given a subgroup $\mathcal{S}$ of the Pauli group, one can define a code from it by considering all the states that are fixed by each element of the group. This indeed constitutes a linear subspace, but in order for it to be nontrivial, one has to restrict $\mathcal{S}$ to be abelian, or commutative, and also require that $-I \notin \mathcal{S}$. In this sense, any abelian subgroup of $\mathcal{P}_N$ that does not contain $-I$ is called a stabilizer group.

One of the advantages of the formalism is that it allows for a compact representation of the code, i.e., it suffices to provide a set of *generators* for the stabilizer group that defines the code. From these *generators*, one can find the codestates (or logical states), in principle, by co-diagonalizing the *generators* and finding their common eigenstates corresponding to the eigenvalue 1. The number $M$ of *generators* necessary to describe the code is related to its dimension by the relation $M = N - K$, much in the way that the number of parity-check equations is related to the code dimension in classical linear codes.

Concerning the encoding procedure, figure 4 presents the example of the five-qubit code, described by *generators* of its stabilizer group, along with the encoding circuit. In the circuit, each wire describes the sequence of gates applied to a qubit. $H$ is the Hadamard gate, given by $H = \frac{1}{\sqrt{2}}Z + \frac{1}{\sqrt{2}}X$, or $|+\rangle\langle 0| + |-\rangle\langle 1|$, which basically swaps the bases of eigenstates of $X$ and $Z$, and is a key ingredient to produce entangled states. Another key ingredient is the
# 1. QUANTUM ERROR CORRECTION

controlled-NOT (CNOT) gate, or controlled-X, denoted by the ⊕ symbol, which acts on a pair of qubits (a control and a target) as $|c\rangle \otimes |t\rangle \mapsto (|0\rangle \langle 0| \otimes I + |1\rangle \langle 1| \otimes X)|c\rangle \otimes |t\rangle$. Other controlled gates can be built, replacing X by other operators, such as the controlled-Y and Z gates that appear in the circuit.

$$ H = \begin{pmatrix} X & Z & Z & X & I \\ I & X & Z & Z & X \\ X & I & X & Z & Z \\ Z & X & I & X & Z \end{pmatrix} $$
$$ \begin{cases} S_1 = X \otimes Z \otimes Z \otimes X \otimes I \\ S_2 = I \otimes X \otimes Z \otimes Z \otimes X \\ S_3 = X \otimes I \otimes X \otimes Z \otimes Z \\ S_4 = Z \otimes X \otimes I \otimes X \otimes Z \end{cases} $$

(a) Stabilizer matrix of the five-qubit code. Each row represents a generator of the stabilizer group, columns represent physical qubits.

<figure><img src="image_6.png" /></figure>

(b) Quantum circuit for encoding the five-qubit code. The lowermost qubit defines the resulting logical state, the remaining are ancillary. Details in [30].

Figure 4: Quantum circuit for encoding the five-qubit code. The lowermost qubit defines the resulting logical state, the remaining are ancillary.

On the decoding side, stabilizer codes share some similarities with their classical counterparts, the linear codes. The stabilizer matrix from figure 4a can be used to measure syndromes just like parity-check equations. More precisely, one can build a quantum circuit that takes as input an N-qubit system prepared in a state $|\psi\rangle$, belonging to the code Q with stabilizer group $S = \langle S_1, \dots, S_M \rangle$. The state is affected by some error $E_\zeta$, which we can assume to be an element of the Pauli group $\mathcal{P}_N$, yielding the modified state $|\tilde{\psi}\rangle = E|\psi\rangle$. The circuit then applies a sequence of gates and then measures, for each generator $S_r$, $r = 1, \dots, M$, the state of an ancilla qubit, obtaining an eigenvalue $\pm 1$, which indicates whether E commutes or anti-commutes with $S_r$. By mapping the eigenvalues to bits as $+1 \mapsto 0, -1 \mapsto 1$, this circuit extracts M classical bits, or $\sigma \in \mathbb{F}_2^M$, which we call the syndrome.

In figure 5 we have a scheme of the syndrome extraction circuit we just described. The uppermost wires are M ancilla qubits, which are prepared in the $|0\rangle$, and then entangled with the other qubits via conjugation by the Hadamard gate and application of controlled gates. The N “data” qubits composing the system are represented by their density operator $\rho$, which undergoes a series of quantum gates, controlled by the ancilla qubits, which are then measured after conjugation by the Hadamard gate. The measurements are with respect to the computational basis, and ultimately yield a syndrome vector $\sigma \in \mathbb{F}_2^M$. The quantum operation described by this circuit can be written as a CPTP map $\mathcal{M}$ given by

$$ \mathcal{M}: \rho \mapsto \sum_{\sigma \in \mathbb{F}_2^M} |\sigma\rangle \langle\sigma| \otimes (\Pi_\sigma \rho \Pi_\sigma^\dagger), \quad \text{with} \quad \Pi_\sigma = \prod_{r=1}^M \frac{I + (-1)^{\sigma_r} S_r}{2}. \qquad (2) $$
<figure><img src="image_7.png" /><figcaption>Figure 5: Syndrome extraction circuit on an N-qubit system, with M ancilla qubits, one for each generator of the stabilizer group S.</figcaption></figure>

Notice that the $\Pi_\sigma$ are the products of the projectors associated to Pauli measurements of each generator $S_r$. The result of this operation is stochastic in general, but when the error operator $E$ is a Pauli operator, the result is deterministic and the syndrome is given by

$$ \sigma_r = \begin{cases} 0, & \text{if } E \text{ and } S_r \text{ commute,} \\ 1, & \text{if } E \text{ and } S_r \text{ anti-commute.} \end{cases} $$

If $E$ is a general unitary error, it can be decomposed as $E = \sum_{P \in \mathcal{P}_N} \alpha_P P$, and then the application of the syndrome extraction circuit to the state $\rho = E |\psi\rangle \langle\psi| E^\dagger$ gives

$$ \mathcal{M}(\rho) = \sum_{\sigma \in \mathbb{F}_2^M} \beta_{\sigma} |\sigma\rangle \langle\sigma| \otimes (P_{\sigma,1} |\psi\rangle \langle\psi| P_{\sigma,2}^{\dagger}), $$

where the $\beta_\sigma$ are obtained from the $\alpha_P$ coefficients, and $P_{\sigma,1}, P_{\sigma,2}$ are elements of the Pauli group that yield syndromes equal to $\sigma$. In summary, the result is a mixed state, but whatever the obtained syndrome is, it determines the resulting state of the system, in the same way it would if the error $E$ was a Pauli operator.

The next step is to decode this syndrome, and this is done classically. Given the syndrome, one needs to infer the error $E$ that occurred. In fact, what is necessary is only a coset $E \cdot S$ of the error, because any two errors that differ by an element of the stabilizer group will have the same effect on any codestate. Errors that belong to the same coset of $S$ are said to be equivalent. One natural strategy to infer the error is to consider the minimum weight decoder, given by

$$ \hat{E} = \arg\min\{P \in \mathcal{P}_N : \sigma(P) = \sigma\}, \quad \text{where} \quad |P| := \#\{j \in \{1, \dots, N\} : P_j \neq I\}. $$

After the syndrome is decoded, the correction is done by applying a Pauli operator $\hat{E}$ that has the same syndrome $\sigma$ as $E$ to the data qubits $\rho$.

Equivalent errors yield the same syndrome, but the converse is not true. In fact, non equivalent errors may yield the same syndrome, just as in the classical setting. In particular, errors that have zero syndrome but are not stabilizers are called logical operators. The set of
# 1. QUANTUM ERROR CORRECTION

logical operators is given by $\mathcal{N}(S) = \{L \in \mathcal{P}_N : LS = SL \ \forall S \in S\}$, i.e., the centralizer of the stabilizer group. Indeed, operators with zero syndrome must commute with all generators of $S$, and so they commute with all of $S$. Their name is due to the fact that these operators act on the code states in $Q$ by permuting them, but leaving the code as whole unchanged. To see this, take any $|\psi\rangle \in Q$ and $L \in \mathcal{N}(S)$, we have

$$L |\psi\rangle = LS |\psi\rangle = SL |\psi\rangle,$$

so $L |\psi\rangle$ is fixed by $S$ for any $S \in S$, thus $L |\psi\rangle \in Q$. The fact that the logical operators are able to alter the logical state of the system makes it clear that $\mathcal{N}(S)$ is strictly larger than $S$. This would only be false for a code of dimension 0. Also, on a side note, it can be shown that $\mathcal{N}(S)$ also coincides with the normalizer of the group $S$, hence the notation.

When they occur as errors, the logical operators are problematic, because they go unnoticed by the syndrome extraction and thus the correction in this case is hopeless. So, it is desirable to have codes for which logical operators are unlikely to occur as errors. Just as in the classical setting, we introduce the notion of distance of a quantum code, which is given by the minimum distance between codestates, or equivalently, the minimum weight of a logical operator, defined as

$$D = \min\{|L| : L \in \mathcal{N}(S) \setminus S\}.$$

Quantum codes are then classified by their parameters as $[[N, K, D]]$, or just $[[N, K]]$ when the distance is unknown.

## 1.3.2 CSS codes

A particular construction of QEC codes was introduced by Calderbank, Shor [31] and Steane [32], and became known as the CSS construction. From the perspective of the stabilizer formalism, CSS codes are simply stabilizer codes for which the stabilizer group is generated by elements of the form $\{\pm 1, \pm i\} \cdot \{I, X\}^{\otimes N}$ or of the form $\{\pm 1, \pm i\} \cdot \{I, Z\}^{\otimes N}$. Notice that the 1-qubit Pauli group is generated by two elements, so we pick $X$ and $Z$ as a sort of basis for it. In terms of syndrome extraction and error correction, everything works in the same way as for general stabilizer codes, the key observation being that this separated structure allows for the independent correction of $X$ errors, known as bit-flips, and $Z$ errors, known as phase-flips. In particular, the $X$-type generators are employed to measure the syndrome of $Z$-type errors, and vice-versa.

General stabilizer codes are represented by their stabilizer matrix as exemplified in figure 4a, which as we have seen also work as parity-check matrices. These matrices should be understood as elements of $\mathbb{F}_4^{M \times N}$, because the additive structure of $\mathbb{F}_4$ is isomorphic to that of the phaseless Pauli group $\tilde{\mathcal{P}}$, or the Klein group. This representation allows for the computation of ranks of such matrices, which are connected to the dimension of the codes. An alternative representation for such matrices is the symplectic representation, given by the mapping

$$ \begin{array}{l} (\cdot)_{\text{sp}} : \tilde{\mathcal{P}}_N \cong \mathbb{F}_4^N \rightarrow \mathbb{F}_2^N \times \mathbb{F}_2^N \cong \mathbb{F}_2^{2N} \\ \bigotimes_{j=1}^N X^{a_j} Z^{b_j} \mapsto ((a_j)_j|(b_j)_j), \end{array} $$
which allows for rank computations in simple $\mathbb{F}_2$ arithmetic. As an example, for the following parity-check matrix its symplectic representation would be

$$H = \begin{pmatrix} X & Z & I & I & Z \\ Z & X & X & Y & I \\ Y & Y & X & Y & Z \end{pmatrix} \iff H_{\text{sp}} = \begin{pmatrix} 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 1 & 1 & 0 \\ 1 & 1 & 1 & 1 & 0 \end{pmatrix} \begin{vmatrix} 0 & 1 & 0 & 0 & 1 \\ 1 & 0 & 0 & 1 & 0 \\ 1 & 1 & 0 & 1 & 1 \end{vmatrix}.$$

In the case of CSS codes formed by $M_X$ X-type generators and $M_Z$ Z-type generators, this representation reduces to the following block structure:

$$H = \left( \begin{array}{c|c} H_X & 0 \\ \hline 0 & H_Z \end{array} \right), \quad \text{where} \quad H_X \in \mathbb{F}_2^{M_X \times N}, H_Z \in \mathbb{F}_2^{M_Z \times N},$$

and for this reason, we can conceive CSS codes as a combination of two classical linear codes, with parity-check matrices $H_X$ and $H_Z$. In order for a pair of classical codes $C_X$ and $C_Z$ to characterize a CSS code, one must require that $C_Z^\perp \subset C_X$, equivalently $C_X^\perp \subset C_Z$. This is because the overall stabilizer group must be abelian, so the X-type and Z-type generators must commute, which can be formulated as $H_X H_Z^T = 0$, equivalently $H_Z H_X^T = 0$.

From the parameters of the classical codes $H_X$, $H_Z$ making up a CSS code, we can compute the quantum code's parameters. If $C_X$, $C_Z$ are $[N, k_X]$, $[N, k_Z]$ codes respectively, the resulting CSS code is a $[[N, K]]$ code with $K = N - \text{rank } H_X - \text{rank } H_Z = k_X + k_Z - N$. As for the distance of the code, we must consider the relation between $C_X$ and $C_Z$ and define

$$D_X := \min\{|e_X| : e_X \in C_X \setminus C_Z^\perp\} \quad \text{and} \quad D_Z := \min\{|e_Z| : e_Z \in C_Z \setminus C_X^\perp\},$$

which yield the minimum distance of the CSS code $D = \min(D_X, D_Z)$.

Just as in the classical case, CSS codes also admit a graphical representation by their quantum Tanner graphs. In figure 6 we have an example of the quantum Tanner graph of the Steane code [32], which is a CSS $[[7, 1, 3]]$ code, produced by taking both $C_X$ and $C_Z$ to be the $[7, 4, 3]$ Hamming code, from figure 2a.

### 1.3.3 The quantum erasure channel

So far we have considered quantum codes for correcting Pauli errors, even if more general errors can also be corrected by the same approach. The most usual noise model for quantum channels is that of the *depolarizing channel*, where each qubit is affected by a Pauli error independently, with probability $p$, split evenly among X, Y and Z type errors:

$$\mathcal{E} : \rho \mapsto (1-p)\rho + \frac{p}{3}X\rho X + \frac{p}{3}Y\rho Y + \frac{p}{3}Z\rho Z.$$

Biased noise models are also conceivable and can be considered realistic depending on the physical platform.

In some platforms, such as Rydberg atom arrays, other kinds of errors can arise such as *leakages*. Such errors may be converted to *erasures* [33], which can be advantageous for the task of error correction.
# 1. QUANTUM ERROR CORRECTION

<figure><img src="image_8.png" /><figcaption>Figure 6: Quantum Tanner graph for the Steane code.</figcaption></figure>

The quantum erasure channel, in analogy to its classical counterpart, is a noise model consisting of an error that deranges an affected qubit in such a way that it becomes manifestly corrupted, and the information about its state is lost. One way to represent that mathematically is to consider that an erased qubit becomes a *maximally mixed state*, that is, a statistical mixture of both basis states in equal proportion, $\rho = \frac{I}{2}$. Alternatively, one can augment the Hilbert space describing a qubit to include a wildcard orthogonal state, modeling it as $\text{span}\{|0\rangle, |1\rangle\} \oplus^{\perp} \mathbb{C}|?\rangle \cong \mathbb{C}^3$. In any case, the premise of the quantum erasure channel is the same as in the classical setting: we are given a random erasure vector $\epsilon \in \mathbb{F}_2^N$ (which is known to the decoder) describing which qubits were erased, while all the remaining qubits are left untouched.

Just as in the classical case, the quantum erasure channel represents a simplified model that facilitates the task of error correction, but also the analysis of quantum codes and decoding algorithms, hence its interest. In order to decode an erasure, the decoder relies on the information provided by the unaffected qubits. If those are enough to determine the codestate unambiguously, the erasure can be corrected.

For an erasure to be non-correctable, there must be an ambiguity over the codestate. In other words, the unaffected qubits must leave room for at least two possible codestates, which are related by some logical operator $L \in \mathcal{N}$. The distinction between the possible codestates must lie in the erased qubits, hence, the support of such logical operator must be contained in the support of the erasure. In the case of CSS codes, a characterization of the correctable erasures is given in [34], Proposition 3.1, which is the pillar of the procedures developed in this work to measure the quality of codes.
**Proposition 1.2.** Let $(C_X, C_Z)$ be classical linear codes defining a CSS code $Q$. An erasure vector $\epsilon \in \mathbb{F}_2^N$ is non-correctable if and only if there is $v \in \mathbb{F}_2^N$ with $\text{supp} v \subset \text{supp} \epsilon$ such that

$$v \in (C_Z \setminus C_X^{\perp}) \cup (C_Z \setminus C_X^{\perp}). \qquad (3)$$

From proposition 1.2 it is possible to derive an algorithm to decide whether an erasure vector is correctable or not. This algorithm will be presented in section 2.

## 1.4 qLDPC and HGP codes

Quantum low-density parity-check codes are defined in a similar way to their classical counterparts, i.e., they are quantum codes whose parity-check matrices have column weights and row weights bounded by constants. In the CSS case, we can define qLDPC codes as codes whose Tanner graphs have bounded degrees. In [35] there is a rich overview of many known families of qLDPC codes, their properties, limitations and details about encoding and decoding procedures. We point out their main property of interest, which is the fact that, because the number of syndrome measurements involving each qubit is bounded, they allow for syndrome extraction circuits of constant depth, which is a key ingredient for achieving fault-tolerant quantum computing with constant overhead [12].

Many constructions for qLDPC codes are known and have interesting properties. In this work, we focus on hypergraph product codes (HGP) [13], which are a way to obtain quantum CSS codes from any two classical codes. The definition of HGP codes relies on the notion of the cartesian product of two graphs. Recall that for two graphs $G_1 = (V_1, E_1)$, $G_2 = (V_2, E_2)$, their cartesian product is given by $G_1 \times G_2 := (V_1 \times V_2, E_\times)$, where

$$E_\times = \{((u, v), (u, v')) : u \in V_1, (v, v') \in E_2\} \cup \{((u, v), (u', v')) : (u, u') \in E_1, (v, v') \in V_2\}.$$

In other words, the product vertices are pairs of vertices from the two factors, and the product edges are drawn between pairs with at most one different component, and whose projections on this differing component are connected in the corresponding factor graph. This construction can be generalized for an arbitrary number of factors.

To build HGP codes, we consider the Tanner graphs of two classical codes $C_1, C_2$, written as $\mathcal{T}_1 = (C_1 \cup B_1, E_1)$, $\mathcal{T}_2 = (C_2 \cup B_2, E_2)$, where $C_i, B_i$ stand for the sets of check nodes and bit nodes, which are disjoint. The HGP code HGP$(C_1, C_2)$ is defined as the CSS code whose quantum Tanner graph is given by $\mathcal{T}_1 \times \mathcal{T}_2$, with its vertex sets $Q \cup C_X \cup C_Z$ organized as

$$\begin{cases} Q = (C_1 \times C_2) \cup (B_1 \times B_2) \\ C_X = C_1 \times B_2, \quad C_Z = B_1 \times C_2 \end{cases}.$$

Notice that since $\mathcal{T}_1$ and $\mathcal{T}_2$ are both bipartite, the resulting graph is multipartite with 4 parts, namely $C_1 \times C_2, C_1 \times B_2, B_1 \times C_2$ and $B_1 \times B_2$, which only have edges among each other, but not internally. The check-check nodes and the bit-bit nodes are taken together to represent the qubits in the HGP code, whereas the other two parts each take the role of either X-type or Z-type quantum checks. Which is which, is up to convention. Figure 7 schematizes this
# 1. QUANTUM ERROR CORRECTION

<figure><img src="image_9.png" /><figcaption>Figure 7: Quantum Tanner graph of the Hypergraph Product code of the [7, 4, 3] Hamming code with itself. Blue and red nodes represent the X and Z check nodes, respectively, while purple nodes represent the qubit nodes.</figcaption></figure>
partition in the case of the code given by the hypergraph product between the [7, 4, 3] Hamming code with itself.

In terms of parity-check matrices, let $H_1 \in \mathbb{F}_2^{m_1 \times n_1}$, $H_2 \in \mathbb{F}_2^{m_2 \times n_2}$ be the parity-check matrices for codes $C_1, C_2$, with parameters $[n_1, k_1], [n_2, k_2]$ respectively, then the CSS parity-check matrices $H_X, H_Z$ of the code HGP($C_1, C_2$) can be written as

$$ \left\{ \begin{array}{l} H_X = \left( \begin{array}{c|c} I_{m_1} \otimes H_2^T & H_1 \otimes I_{n_2} \\ \hline H_1^T \otimes I_{m_2} & I_{n_1} \otimes H_2 \end{array} \right) \in \mathbb{F}_2^{M_X \times N}, \\ H_Z = \left( \begin{array}{c|c} H_1^T \otimes I_{m_2} & I_{n_1} \otimes H_2 \\ \hline H_2^T \otimes I_{n_1} & I_{n_2} \otimes H_1 \end{array} \right) \in \mathbb{F}_2^{M_Z \times N}, \end{array} \right. \quad (4) $$

where $M_X = m_1n_2$, $M_Z = n_1m_2$, and $N = m_1m_2 + n_1n_2$. Depending on the convention and the ordering of the qubits, the formulas in (4) may change, as in [13] and [26]. Nevertheless, the definition remains equivalent. One can verify that indeed, $H_X$ and $H_Z$ define a CSS code by the computation

$$ H_X H_Z^T = (I_{m_1} \otimes H_2^T \quad H_1 \otimes I_{n_2}) \begin{pmatrix} H_1 \otimes I_{m_2} \\ I_{n_1} \otimes H_2^T \end{pmatrix} = H_1 \otimes H_2^T + H_1 \otimes H_2^T = 0, $$

since we are working in $\mathbb{F}_2$.

Already we can see how the number of physical qubits $N$ relates to the parameters of the two classical codes, and further relations exist for the other parameters. In summary, the HGP of two classical codes with parameters $[n_1, k_1, d_1]$, $[n_2, k_2, d_2]$, assuming they are represented by full-rank parity-check matrices, gives a CSS code with parameters $[[m_1m_2 + n_1n_2, k_1k_2, d_Q]]$, where $d_Q \ge \min(d_1, d_2)$. In general, the HGP construction depends on the particular parity-check (or equivalently, the Tanner graph) representation of the code. If the matrices are not full-rank, the dimension of the resulting HGP code is given by $k_1k_2 + k_1^T k_2^T$, where $k_i^T$ refers to the dimension of the linear code defined by $H_i^T$. When $H_i$ is full-rank, $k_i^T = 0$, so this expression generalizes the formula for the dimension of the HGP code. As for the distance of the code, the lower bound can be an equality in some cases where one of $H_i$ is not full-rank. Moreover, the LDPC property is also preserved by the HGP construction, i.e., if the two classical codes $C_1$ and $C_2$ are LDPC, so is the quantum code HGP($C_1, C_2$). The relation between the bounds on the degrees can be deduced by looking at the degrees in the classical Tanner graphs and considering the connectivity structure of their cartesian product.

By considering HGP products of a single classical code with itself, one can prove the existence of families of HGP codes with the asymptotic parameters $[[N, \Theta(N), \Theta(\sqrt{N})]]$, which is optimal in terms of dimension, i.e., the family has constant rate $K/N$. In terms of minimum distance, better codes have been found [15, 16], achieving both $K = \Theta(N)$ and $D = \Theta(N)$. For a summary of those codes, visit [35]. The recent work [36] carries out a comparison of the performance of many kinds of product codes, i.e., CSS codes produced from classical codes via some sort of product construction.
# 2 Optimizing QEC codes

All the code containing the implementations and results of the experiments is available the project's github repository.

In this work, we are focused on optimizing HGP codes. By adopting this particular construction, the problem of optimizing the QEC code is reduced to that of optimizing the two classical codes making up the hypergraph product. In particular, we consider only HGP codes of a single classical code $C$ with itself, of the form $HGP(C, C)$, or $HPG(C)$ for short. This restriction induces a symmetry between the $C_X$ and $C_Z$ codes making up the quantum code, i.e., they are isomorphic. This is a reasonable property when considering noise models where bit-flips and phase-flips are equally likely. For biased noise models, it can be interesting to consider asymmetric codes, built from two different classical codes. Even though the approach developed here can be adapted to this asymmetric case, the consideration of only HGP codes built from a single classical code with itself simplifies the modeling of the optimization problem, and by consequence, greatly reduces the computational effort of the optimization methods.

One motivation for the task of optimizing such codes comes from [4]. In that work, the authors develop a new decoder for HGP codes under the quantum erasure channel, the VH decoder. The objective of such decoder is to be efficient, while keeping a performance comparable to that of the optimal, maximum likelihood decoder, based on Gaussian elimination. The analysis of the weak points of this VH decoder indicates that its stopping sets, that is, sets of bits such that, when erased, lead to decoding failure, are composed of stopping sets in the classical codes making up the hypergraph product. In this context, optimizing the classical codes to avoid structures that produce those stopping sets might be the key to producing codes that perform better, not only for the optimal decoder, but also for the faster, VH decoder.

In [4], the method used to generate the classical codes making up the HGP codes is the progressive edge growth algorithm [17], which generates Tanner graphs of given dimensions and degree distributions (i.e., produce an LDPC code), by trying to maximize their girth, that is, the length of the shortest cycle in the graph. Heuristically, short cycles in the Tanner graph lead to bad performance, so it makes sense to adopt such strategy. However, it remains an ad hoc method, as the precise characteristics of what makes a qLDPC code optimal are not yet well understood.

### 2.1 Modeling the problem in RL framework

As an alternative to classical ad hoc methods, this work adopts a search-based approach to find better codes among the class of HGP codes.

Inspired by [23], we frame our problem as that of searching for a good quantum code inside a collection of quantum codes, which will further be referred to as a *state space*, as in the jargon of Reinforcement Learning frameworks. This search is realized through the application of *actions*, that consist on *locally* modifying the codes, thus obtaining new ones, effectively traversing the state space through such actions. In order to find better codes, one has to define a *cost function* with which to measure the quality of such codes. Finally, an algorithm is needed to optimize the codes under such cost function over the state space.
In [23], the authors consider a state space comprised of surface codes [24], in which the actions consist of local deformations, by adding or removing physical qubits to the patch defining the surface code. Their cost function is based on a benchmarking tool called SQUAB [25], which runs in linear time (in terms of the number of physical qubits $N$) and assesses the performance of the surface code under an optimal decoder for the quantum erasure channel. They then apply the Projective Simulation (PS) algorithm [21, 22] to optimize the codes.

In our context, many adaptations are required. To fix notation, we'll denote the state space by $\mathcal{S}$, and for a state $s \in \mathcal{S}$, its action space will be denoted $\mathcal{A}_s$. Actions $a \in \mathcal{A}_s$ are directly connected with state *transitions*. As opposed to general Markov Decision Processes (MDPs) used as a framework for RL, here the transitions are deterministic given the action. In other words, to each action $a_{s'} \in \mathcal{A}_s$, is associated a single *next state* $s' \in \mathcal{S}$. In this context, whenever we talk about *reward functions*, it's equivalent to associate them to either a state-action pair $(s, a) \in \mathcal{S} \times \mathcal{A}_s$ or to a state-transition pair $(s, s') \in \mathcal{S} \times \mathcal{S}$, where $s'$ can be reached from $s$ by some action $a \in \mathcal{A}_s$.

### 2.1.1 State space

For our setting, we want to consider a state space comprised by classical LDPC codes with $m$ checks and $n$ bits, and a given pair of left and right-degree distributions, $(l_i)_{i=1}^m$, $(r_j)_{j=1}^n$. The degree distributions must of course be compatible, since by double-counting the edges of the Tanner graph (or the 1s in the parity-check matrix) one must have $\sum_i l_i = \sum_j r_j = |E| = \sum_{i,j} H_{ij}$, where $\mathcal{T}(\mathcal{C}) = (C \cup B, E)$ is the Tanner graph of the code, and $H \in \mathbb{F}_2^{m \times n}$ is its parity-check matrix.

In particular, we are interested in *regular* or almost regular codes, i.e., codes for which the degree distributions $(l_i)_{i}$, $(r_j)_{j}$ are close to uniform. There is a systematic way to generate regular LDPC codes by first picking the degrees $l, r$. The number of edges $E$ will then be a multiple of $lcm(l, r)$, say $E = k \cdot lcm(l, r)$. Next, build a bipartite graph with two vertex sets of size $E$, choose a permutation of $E$ elements and connect each left node to the right node corresponding to the permutation. Finally, pose $m = E/l, n = E/r$, and group the left vertices by clusters of size $l$, and the right vertices by clusters of size $r$. This will yield a bipartite graph with $m$ left and $n$ right nodes, with respective left and right degrees equal to $l$ and $r$. This construction however is prone to yield graphs with parallel edges, or *multigraphs*.

For the sake of technical simplicity, we will model the state space $\mathcal{S}$ as the set of bipartite multigraphs with vertex sets of fixed sizes $m, n$ with fixed left and right degree distributions $(l_i)_{i}$, $(r_j)_{j}$. However, when it comes to evaluating the cost function on a given state, we will consider the underlying *simple graph*, that is, we take the multigraph and simply discard any repeated edges. This will produce a simple Tanner graph that properly defines a linear code. Notice, though, that many different graphs can be *isomorphic*, meaning, they are equal up to a permutation of the vertices. Isomorphic graphs induce the exact same codes, up to reordering the bits. Moreover, from the generator matrix point of view, replacing any row by a linear combination of the others still describes the same code, even if the Tanner graphs are not isomorphic. For technical reasons, we did not adopt any form of equivalence identification, though the idea was investigated. More details in 3.4.
## 2. OPTIMIZING QEC CODES

Another reason for why it is simpler to model the state space like that will become clear when we discuss the action space.

### 2.1.2 Action space

The choice of our action space is arguably the most important part of the modeling. The state space must be closed under the actions we define, and ideally it should also be connected. The main motivation behind our choice, though, is the expected effect the actions may have on the performance of a code.

For each state $s \in S$, the action space $A_s$ will be given by the set of pairs of edges from the bipartite multigraph $s = (A \cup B, E)$. The concrete action consists of taking any two edges $(u_1, v_1), (u_2, v_2) \in E$, and replacing them by the crossed edges $(u_1, v_2), (u_2, v_1)$. In terms of the states, if $s = (A \cup B, E)$, the new state is given by $s' = (A \cup B, E')$, where

$$E' = (E \setminus \{(u_1, v_1), (u_2, v_2)\}) \cup \{(u_1, v_2), (u_2, v_1)\}.$$

<figure><img src="image_10.png" /><figcaption>Figure 8: Illustration of the edge-crossing operation. The plots represent states $s = G_1, s' = G_2, s'' = G_3$, and the highlighted edges are the ones being crossed at each transition.</figcaption></figure>

Figure 8 illustrates this operation. The state $s = G_1$ is the Tanner graph of the [7, 4, 3] Hamming code, and upon application of the edge-crossing action to the edges $(c_1, n_5), (c_2, n_6)$,
the resulting state $s' = G_2$ has them replaced by the new edges $(c_1, n_6)$, $(c_2, n_5)$. On a second application of the operation, the state $s'$ is modified to obtain $s'' = G_3$, this time applying the edge-crossing to another pair of edges, $(c_1, n_2)$, $(c_3, n_5)$.

Notice that because we allow for parallel edges, $E$ and $E'$ should be understood as multisets, i.e., they may hold copies of a same edge, and the replacement operation acts only on one copy the involved edges. Also notice that in the case where $u_1 = u_2$ or $v_1 = v_2$, the edge-crossing operations doesn't change the graph, i.e., $s' = s$.

To implement such action, we enumerate the edges of $s$, $E = \{e_j\}_{j=1}^{|E|}$, sorted lexicographically, and then index the actions by a pair of indices $1 \le i < j \le |E|$, indicating the edges to be crossed. This pair can be converted into a single index via the formula

$$a = j - i + \sum_{k=1}^{i-1} |E| - k = j + (i-1)|E| - \binom{i+1}{2} \in \{1, \dots, \binom{E}{2}\},$$

and conversely

$$\left\{ \begin{array}{l} i = \left\lfloor \frac{2|E| + 1 - \sqrt{(2|E| - 1)^2 - 8(a - 1)}}{2} \right\rfloor \\ j = a - (i - 1)|E| + \binom{i+1}{2} \end{array} \right. .$$

### 2.1.3 Cost function

In order to measure the performance of each state (or quantum code), we chose a similar approach as that of [23], that is, by doing a Monte Carlo estimation of the logical error rate $p_L$ of the code under the quantum erasure channel.

Algorithm 1 computes the plain Monte Carlo estimation of the logical error rate based on the correctability criterion from proposition 1.2. For CSS codes, as is the case of HGP codes, the criterion to test whether an erasure vector $\epsilon \in \mathbb{F}_2^N$ is correctable reduces to a linear algebraic problem.

Let $H_X$, $H_Z$ be the parity-check matrices of the quantum code, and for any matrix $A \in \mathbb{F}_2^{M \times N}$ denote by $A^\epsilon$ the matrix obtained by extracting the columns of $A$, $(A_{.,j})_j$ where $\epsilon_j = 1$. In other words, select the columns associated to the erased qubits. Analogously, denote by $A^{\bar{\epsilon}}$ the matrix obtained from the non-erased columns. For vectors $v \in \mathbb{F}_2^N$, seen as column matrices, denote by $v_\epsilon$ the extraction of the erased rows, and $v_{\bar{\epsilon}}$ the non-erased rows.

As in proposition 1.2, $\epsilon$ is non-correctable when $\exists v \in \mathbb{F}_2^N$ with $\text{supp } v \subset \text{supp } \epsilon$ satisfying $v \in \text{ker } H_Z \setminus \text{im } H_X^T$ or $v \in \text{ker } H_X \setminus \text{im } H_Z^T$. The condition $\text{supp } v \subset \text{supp } \epsilon$ implies that we can parameterize $v$ by $v_\epsilon$, because $v_{\bar{\epsilon}} = 0$. Up to reordering the qubits, we have

$$v = \binom{v_{\epsilon}}{0}, \quad H_X = \binom{H_X^{\epsilon}}{H_X^{\bar{\epsilon}}}, \quad H_Z = \binom{H_Z^{\epsilon}}{H_Z^{\bar{\epsilon}}}. $$

Now the condition $v \in \text{ker } H_Z$ reduces to $v_\epsilon \in \text{ker } H_Z^\epsilon$. In order to test the condition $v \notin \text{im } H_X^T$ given this parametrization, it is convenient to consider parity-check matrix for which $H_X^T$ is the corresponding generator matrix. In the notation of lemma 1.1, we have

$$v \notin \text{im } H_X^T \iff v \notin \text{ker } \eta(H_X^T),$$
## 2. OPTIMIZING QEC CODES

and by adopting the previous parametrization, we have

$$v_{\epsilon} \notin \ker \eta(H_{X}^{T})^{\epsilon}.$$

Now, in order to check for the existence of $v_{\epsilon} \in \ker H_{Z}^{\epsilon}$ such that $v_{\epsilon} \notin \ker \eta(H_{X}^{T})^{\epsilon}$, we have to check the condition for all vectors in a basis of $\ker H_{Z}^{\epsilon}$. For this we consider a generator matrix for $\ker H_{Z}^{\epsilon}$, which we denote by $\gamma(H_{Z}^{\epsilon})$. With the established notation, the correctability condition for $\epsilon$ translates to:

$$ \epsilon \text{ is correctable} \iff \left\{ \begin{array}{l} \eta(H_X^T)^\epsilon \gamma(H_Z^\epsilon) = 0 \\ \eta(H_Z^T)^\epsilon \gamma(H_X^\epsilon) = 0 \end{array} \right. . \qquad (5) $$

**Algorithm 1** Monte Carlo estimation of the failure rate for the erasure channel.

**Input:** Number of trials $M$, physical erasure rate $p$, CSS matrices $H_X, H_Z$.

**Output:** Estimated logical error rate $\hat{p}_L$.

1: **Initialization:**

2: Pre-compute matrices $\eta(H_X^T), \eta(H_Z^T)$.

3: $n \leftarrow 0$

4: **for** $i \leftarrow 1, \dots, M$ **do**

5: &nbsp;&nbsp;&nbsp;&nbsp;Sample an erasure vector $\epsilon \sim \mathcal{B}(p)^{\otimes N}$.

6: &nbsp;&nbsp;&nbsp;&nbsp;Compute matrix products $\eta(H_X^T)^\epsilon \gamma(H_Z^\epsilon)$ and $\eta(H_Z^T)^\epsilon \gamma(H_X^\epsilon)$.

7: &nbsp;&nbsp;&nbsp;&nbsp;**if** $\eta(H_X^T)^\epsilon \gamma(H_Z^\epsilon) \neq 0$ or $\eta(H_Z^T)^\epsilon \gamma(H_X^\epsilon) \neq 0$ **then**

8: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$n \leftarrow n + 1$

9: **end if**

10: **end for**

11: **return** $\hat{p}_L = \frac{n}{M}$.

Notice that since algorithm 1 does a plain Monte Carlo simulation, the precision with which we can estimate the failure rate is limited by the number of trials. This can prove to be a problem when we want to measure the failure rate of codes with very good performance, or for very low erasure rates. In fact, more advanced simulation techniques could be used to alleviate this constraint, such as *importance sampling*. However, my attempts to make it work were not fruitful, yielding estimators that systematically underestimated the failure rate and whose variance did not seem to be reduced in comparison to plain MC. Throughout all the experiments, the cost function used to evaluate the codes was given by the result of algorithm 1 with a fixed erasure rate, chosen in a way that, for a fixed budget $M$, the estimate of the failure rate had a variance below 5% of the mean value.

The operators $\eta$ and $\gamma$ are implemented relying on the PLUQ factorization in $\mathbb{F}_2$. To compute $\gamma(H)$, it suffices to solve the linear system $Hx = 0$ and assemble all the linearly independent solutions into a matrix. To compute $\eta(G)$, one relies on the observation that, if $G$ is of the form

$$G = \left( \frac{I}{A} \right), \quad \text{then} \quad H = \left( A \mid I \right) \text{ satisfies } HG = A + A = 0.$$
With this observation, one can simply compute the row-reduced echelon form of $G^T$ and re-order its columns to obtain a form $(I \ A^T)$. After that, applying the corresponding column permutation to $(A \ I)$ yields $H = \eta(G)$.

The implementation of all the linear algebraic operations in $\mathbb{F}_2$ were done with the help of the M4RI library, [37]. It is an extremely optimized C library for doing linear algebra with dense matrices over $\mathbb{F}_2$. To integrate it with the Python code of the project, a C extension for Python was developed.

## 2.2 Projective Simulation

The main approach exploited in this work for the task of optimizing QEC codes is based on the Projection Simulation (PS) [21, 22] framework. Python implementations of the framework and some generalizations are available at [38]. We will present here a context-tailored description of the framework, including its main features. As a RL framework, it models the interaction between an *agent* and an *environment*, where the environment is responsible for issuing *percepts* or *observations* of the underlying state, as well as rewards for the actions taken by the agent.

In the modeling of the QEC optimization problem, actions are transitions, and the rewards are functions of the newly reached state. We can thus think of the combination of state and action spaces as a sort of meta graph, where each state is a vertex and the actions are edges. The cost function (and by consequence, the reward function) can be seen as a signal on this meta graph. Our goal is to find the optimum of this signal by traversing the meta graph in a controlled way.

The PS agent's behavior is governed by local stochastic policies, that allow for a certain amount of exploration. Specifically, for each state $s$, the agent holds a probability distribution $\alpha$ over the action space $\mathcal{A}_s$. In other words,

$$ \mathbb{P}[\text{agent takes action } a_j | \text{agent is in state } s_i] = \alpha_{ij}. $$

The learning algorithm uses the rewards obtained along the exploration to reinforce actions that lead to the accumulation of such rewards, by updating these probability distributions for all states at each step of the exploration.

The concrete implementation of the PS agent consists on two matrices, $h, g \in \mathbb{R}^{|S| \times |A|}$. Each row in $h$ holds the local policy of a state, which is converted into a probability distribution via a softmax transformation:

$$ \alpha_{ij} \propto e^{\beta h_{ij}}, \quad \text{with normalization factor} \quad \sum_{k=1}^{|A_{s_i}|} e^{\beta h_{ik}}, \qquad (6) $$

where $\beta$ is a hyper-parameter of the model. In the initialization, the $h$ matrix is set to zero, i.e., all the distributions $\alpha$ are uniform.

The matrix $g$ is called the *glow matrix*, in analogy to a sort of memory the agent has of its previous actions. It is used to update the $h$ matrix upon reception of a reward $\lambda$ by the formula

$$ h \leftarrow h + \lambda \cdot g. $$
## 2. OPTIMIZING QEC CODES

On top of that, the PS agent also possesses two forgetting mechanisms, controlled by parameters $\gamma \in [0, 1]$, called the forgetting parameter and $\eta \in [0, 1]$, called the glow damping parameter. They act on the $h$ and $g$ matrices respectively by shrinking their weights:

$$h \leftarrow (1 - \gamma)h, \quad g \leftarrow (1 - \eta)g.$$

The forgetting update takes place right before the incorporation of a reward, so the full update of the $h$ matrix takes the form

$$h_{ij}^{(t+1)} = (1 - \gamma)h_{ij}^{(t)} + \lambda^{(t)} \cdot g_{ij}^{(t)}. \qquad (7)$$

As for the glow matrix, in addition to the glow damping mechanism, one entry is set to 1 right after an action is taken:

$$g_{ij}^{(t+1)} = \max \left( (1 - \eta)g_{ij}^{(t)}, \mathbf{1}[s^{(t)} = s_i] \mathbf{1}[a^{(t)} = a_j] \right). \qquad (8)$$

The training algorithm of the PS agent is summarized in algorithm 2, describing the outline of a training episode. The full training may comprise several episodes, resetting the state to $s^{(0)}$ but keeping the $h$ and $g$ weights across episodes.

**Algorithm 2** Training episode in Projective Simulation

**Input:** Hyper-parameters $\gamma$, $\eta$ and $\beta$, initial state $s^{(0)}$.

**Output:** Trained policy $h$.

1: **Initialization:**

2: Initialize $h^{(0)} = 0, g^{(0)} = 0$.

3: **Repeat:**

4: Sample a random action $a^{(t)} \sim \alpha_i^{(t)}$, where $\alpha_i^{(t)} \propto e^{\beta h_i^{(t)}}$, for $s_i = s^{(t)}$.

5: Apply action $a^{(t)}$ and observe a new state $s^{(t+1)}$.

6: Environment issues a reward $\lambda^{(t)} \ge 0$.

7: Update weights $h_{ij}^{(t+1)} \leftarrow (1 - \gamma)h_{ij}^{(t)} + \lambda^{(t)} \cdot g_{ij}^{(t)}$.

8: Update glow matrix $g_{ij}^{(t+1)} = \max \left( (1 - \eta)g_{ij}^{(t)}, \mathbf{1}[s^{(t)} = s_i] \mathbf{1}[a^{(t)} = a_j] \right)$.

9: **Until:** End of the episode.

### 2.2.1 Implementation details

To apply Projective Simulation to our problem, some technical aspects need to be addressed. For most of the project's implementation, we made use of the NetworkX [39] library, as well as other more commonplace Python libraries, such as NumPy [40] and SciPy [41].

First of all, the state is stored by the environment as a multigraph. Performing the edge-crossing operation on it is straightforward with the help of the formulas in section 2.1.2.

Considering the potential size of the state space [42], the agent's matrices $h$ and $g$ are stored as sparse matrices. This is done because, for instance, $g$ only has nonzero entries for each pair of state and action that it has visited along the training. Since $h$ is also initialized with zeros,
it only gets nonzero entries when a nonzero reward is received, so both matrices are sparse. Moreover, in order to index the matrices by the corresponding states, a dictionary is kept to allow for an efficient association between a state and its corresponding row index.

Thanks to the design choice of allowing for multigraphs, the number of edges for any state is the same. By allowing actions that eventually yield to the same state, or that create parallel edges, we manage to have action spaces of the same size for each state, which is convenient for the encoding of actions and for the storage of the matrices h and g.

Concerning the action sampling procedure, since we store each row of the h matrix in sparse format, a clever sampling algorithm is necessary.

* **Stratified sampling:** The idea is that we have a random variable $a \sim \mathcal{P}(\{1, \dots, A\})$ representing the index of the sampled action, A being the total number of actions, with distribution given by $\mathbb{P}[a = j] = p_{ij} \propto \exp(\beta h_{ij})$. However, since we want to store only the nonzero entries of $(h_{ij})_j$, we cannot sample it directly from the vector of probabilities $(p_{ij})_j$. Moreover, sampling from a custom vector of probabilities has time (and space) complexity proportional to the vector size, in this case, $O(A)$. Thus, instead of sampling directly from $(p_{ij})_j$, we'll consider a Bernoulli variable $\xi \sim \mathcal{B}(q)$, a uniform variable $u \sim \mathcal{U}(\{1, \dots, A\})$ and a multinomial variable $\tilde{a}$ with the same support as $(h_{ij})_j$, all independent. The idea is that we can choose the parameter $q$ and the distribution of $\tilde{a}$ such that we have $a = \xi \tilde{a} + (1 - \xi)u$.

* The derivation is as follows: let $\mathcal{H} := \text{supp}(h_{ij})_j$, so

$$ \mathbb{P}[a = j] = \mathbb{P}[\xi = 1] \cdot \mathbb{P}[\tilde{a} = j] + \mathbb{P}[\xi = 0] \cdot \mathbb{P}[u = j]. $$

If $j \notin \mathcal{H} = \text{supp} \tilde{a}$, we have $\mathbb{P}[a = j] = (1 - q) \frac{1}{A}$. Now, since we want

$$ \mathbb{P}[a = j] = \frac{1}{A - |\mathcal{H}| + \sum_{k \in \mathcal{H}} \exp(\beta h_{ik})}, $$

we must then take $q = \frac{\sum_{k \in \mathcal{H}} [\exp(\beta h_{ik}) - 1]}{A - |\mathcal{H}| + \sum_{k \in \mathcal{H}} \exp(\beta h_{ik})}$. Finally, for $j \in \mathcal{H}$, we have

$$ \mathbb{P}[a = j] = \frac{\sum_{k \in \mathcal{H}} [\exp(\beta h_{ik}) - 1]}{A - |\mathcal{H}| + \sum_{k \in \mathcal{H}} \exp(\beta h_{ik})} \cdot \mathbb{P}[\tilde{a} = j] + \frac{1}{A - |\mathcal{H}| + \sum_{k \in \mathcal{H}} \exp(\beta h_{ik})} \\ = \frac{\exp(\beta h_{ij})}{A - |\mathcal{H}| + \sum_{k \in \mathcal{H}} \exp(\beta h_{ik})}, $$

so it suffices to take $\mathbb{P}[\tilde{a} = j] \propto \exp(\beta h_{ij}) - 1$.

Finally, concerning the reward function, it can be chosen to be an encapsulation of the underlying cost function. In [23], the cost function is computed by SQUAB [25] as an estimate of the logical error rate of the code, but the reward function is a sparse version of it, giving a 0 reward if the failure rate is above a fixed threshold, and 1 if it's below the threshold. In this work, we apply the same strategy, by adjusting the threshold as a hyper-parameter.
### 2.3 Simulated Annealing and plain exploration

To assess the applicability of Projective Simulation to the QEC optimization problem, we also implemented an alternative optimization method, called *Simulated Annealing* [43]. Simulated Annealing (SA) is an iterative, heuristic algorithm that exploits stochasticity to approximate solutions for optimization problems.

The algorithm consists on randomly exploring the variable space, i.e., performing a random walk, while guided by a temperature scheduling function. At each transition, the algorithm samples a random neighbor point, and compares its cost function with that of the current point. If the neighbor's value is better, the transition is accepted. However, the key ingredient is that when the neighbor's value is worse, the algorithm can randomly accept the new value and realize the transition all the same. The probability of transition to a worse state depends on the temperature as

$$P_{accept} := \exp(-\Delta/T), \qquad (9)$$

where $T$ is the current temperature, and $\Delta$ is difference between the current value and the neighbor's value. The name *annealing* comes from the fact that the temperature is scheduled to some high value at the beginning $T_0 > 0$, and gradually decays as the iterations progress. So, in the beginning, the temperature is high, and thus the algorithm is highly susceptible to accepting points with worse values, but as the iterations progress, it becomes colder, and the algorithm tends to accept only solutions that improve the function value. The heuristics behind SA is to provide some exploratory behavior in the beginning, in hopes of avoiding local minima, and cooling down in the end to provide some stability.

In our problem formulation, SA can be readily applied, by taking the random neighbor generation to be restricted to neighbors in the meta graph, i.e., states that can be reached from the current state by applying the actions defined in 2.1.2. Throughout the experiments, we adopted a temperature scheduling of the form

$$T = \frac{1}{1 + \beta(t/t_{max})^2}, \qquad (10)$$

where $t$ is the iteration number, $t_{max}$ is the maximum number of iterations, and $\beta$ is a hyper-parameter.

Yet another alternative to Projective Simulation is just a simple exploration of the state space, navigated through by the application of actions. We did experiments with a simple exploration algorithm that starts from a given initial state, and performs a random walk of length $L$. At each step, it computes the cost function of the current state, and also that of a number $N$ of its neighbors. The objective of this exploration is to have a notion of how the cost function varies across the state space, and among neighboring states. Due to the sheer number of different states, it is impractical to perform a full exploration of the space when the codes involved are large.
# 3 Results and Discussion

Along all the experiments we carried out, we focused on four code families, which are the ones obtained by applying the actions to four initial codes. These initial codes are the ones used in [4], and are available at the code repository associated to that paper. The repository stores the classical LDPC codes which are used to produce the HGP codes. These codes have different sizes, and are all produced using the progressive edge-growth algorithm [17] to yield almost regular, (3, 4)-LDPC codes. The resulting HGP codes obtained from them have parameters [[625, 25]], [[1225, 65]], [[1600, 64]] and [[2025, 81]].

One important **disclaimer** must be made about the second code family [[1225, 65]]. Notice that the dimension of the initial code is not a perfect square, which means that the classical code used to produce it is not represented by a full-rank parity-check matrix. During the traversal of the state space, this code is deformed and can have its rank augmented, thus lowering the dimension of the resulting HGP code. Hence, whenever the code family [[1225, 65]] is mentioned in the results, bear in mind that the dimension of the codes in the family may be smaller, down to at least 49, and this is expected to be traded off for an improvement in terms of logical error rate.

Because the complexity of estimating the cost function grows with the number of qubits as $O(N^3)$, the number of Monte Carlo trials was adjusted to yield a reasonable computation time for the experiments. For all the optimization experiments, we fixed a Monte Carlo budget of $M = 10^4$ trials, and fixed the erasure rates for the experiments with each code family to the values in table 1. The choice of such values was based on the relative precision of the estimation of the corresponding failure rate with the predefined Monte Carlo budget. We aimed at a ratio standard deviation/mean of at most 5%.

<table><caption>Table 1: Chosen erasure rates for the optimization experiments on each code family.</caption><thead><tr><th>Code family</th><th>[[625, 25]]</th><th>[[1225, 65]]</th><th>[[1600, 64]]</th><th>[[2025, 81]]</th></tr></thead><tbody><tr><th>Erasure rate</th><td>9/32</td><td>8/32</td><td>9/32</td><td>12/32</td></tr></tbody></table>

Figure 9 shows the performance of each of these codes as measured by the condition in proposition 1.2, i.e., the failure rate of the CSS code under the quantum erasure channel, for a range of erasure rates. The failure criterion adopted is stricter than the failure rate of an actual decoder, and it takes into account both *X* errors and *Z* errors, therefore, the exact values in figure 9 will differ from those presented in [4].

## 3.1 Plain Exploration

For the experiments with plain exploration, we chose a maximum number of codes to be evaluated inside each family, starting the exploration from the initial codes. For each family, we fix a number *N* of neighbors to explore at each step of the random walk of fixed length *L*. In table 2 we have the relation of parameters for each code family, as well as the results of the best codes found the exploration.
## 3. RESULTS AND DISCUSSION

<figure><img src="image_11.png" /><figcaption>Figure 9: Estimation of the performance of the CSS codes from [4]. The number of Monte Carlo trials was chosen to balance the computational cost and the precision of the estimates.</figcaption></figure>

<table><caption>Table 2: Parameters for the plain exploration on each family.</caption><thead><tr><th>Code family</th><th>[[625, 25]]</th><th>[[1225, 65]]</th><th>[[1600, 64]]</th><th>[[2025, 81]]</th></tr></thead><tbody><tr><td><i>N</i></td><td>24</td><td>15</td><td>12</td><td>8</td></tr><tr><td><i>L</i></td><td>120</td><td>70</td><td>40</td><td>30</td></tr><tr><td>Best performance gain</td><td>3.19</td><td>60.00</td><td>31.83</td><td>2.16</td></tr></tbody></table>

In figure 10 we plot the estimated logical error rate for each state visited during the exploration, for each code family. Red dots are the neighbor points, which are probed but don't necessarily belong to the path of the random walk. In blue is the curve of performance for the points belonging to the path of the random walk. The green star represents the best code found along the exploration.

<figure><img src="image_12.png" /><figcaption>Figure 10: Random exploration of the cost function along the state space for each code family. The blue line is the random walk path, and red dots are neighbor points. The green star identifies the best code found. Cyan ticks depict the 95% confidence intervals for each measurement.</figcaption></figure>

As the plots show, the variation of the cost function oscillates pretty wildly across the state space, far beyond the variance of the estimates. In particular, we point out that most of the codes, even the best ones, have neighbors of very poor performance. This sheds some light on the effect of the actions defined in 2.1.2 over the codes' performance.

Another aspect of the behavior of the cost function across the state space can be seen in figure 11, where we can see the distribution of the values of the cost function over the
sampled points. Notice that the code families [[1225, 65]] and [[1600, 64]] have plenty of room
for improvement over the initial codes. For families [[625, 25]] and [[2025, 81]] the performance
gains are less impressive. In the case of [[2025, 81]] this could be simply due to insufficient
exploration, however in [[625, 25]] apparently there isn't much room for improvement, as it was
much more explored than the other families.

<figure><img src="image_13.png" /><figcaption>Figure 11: Histogram of the distribution of the cost function values across the explored states. The black line denotes the performance of the initial code in each family.</figcaption></figure>

## 3.2 Simulated Annealing

Based on the results of the exploration experiments, we have a notion of how much the cost
function varies across the state space and also some baselines for the performance gains.

### 3.2.1 Plain acceptance probability function

With Simulated Annealing, we did two sets of experiments. In the first set, we ran the SA
algorithm for each code family with a maximum number of iterations $t_{max}$, using the tempera-
ture scheduling given in (10), and varying the $\beta$ parameter between 5, 10, 25, 50, 100, 200, 400.
In figure 12, we show the history of the values of the states visited by the algorithm. For
compactness, we display only the top 4 $\beta$ values for each code family.

The first important observation we can draw from figure 12 is that Simulated Annealing
doesn't always converge to optimal codes. In reality, in most simulations the best code found
during the algorithm appears in the middle of the execution, and the final code rarely performs
as well as the best code.

After picking the best codes found during the optimization, we compared their performances
with the initial code in each family over a range of erasure rates, just as in figure 9. The plots
in figure 13 show that codes that were optimized for a particular noise level also perform better
in a range of other noise levels. Another noticeable fact is that the intermediate sized families
[[1225, 65]] and [[1600, 64]] exhibit dramatic improvements over the initial codes, with close to
two orders of magnitude of performance gain in some noise levels.
## 3. RESULTS AND DISCUSSION

<figure><img src="image_14.png" /><figcaption>(a) Evolution of the code performances during SA in the code family [[625, 25]].</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(b) Evolution of the code performances during SA in the code family [[1225, 65]].</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(c) Evolution of the code performances during SA in the code family [[1600, 64]].</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(d) Evolution of the code performances during SA in the code family [[2025, 81]].</figcaption></figure>

Figure 12: Evolution of the code performances throughout the execution of Simulated Annealing in each code family, for different values of $\beta$.

### 3.2.2 Log-scale acceptance probability function

The second set of experiments were carried out in the same fashion as the first one, but now we changed the acceptance probability function of (9) to take into account the variation of the logarithm of the cost function, instead of the plain cost function.

Studying the general trend of the curves in both figures 12 and 14, we can see that larger $\beta$ values tend to produce a decreasing trend in the history of the cost function. However, this effect is more pronounced in figure 14, where the acceptance probability function works in log-scale.

In fact, from table 3 we can see that this parametrization also yields better results in terms of best performance gain for each code family. The whole motivation behind choosing the log-scale acceptance function is related to how the cost function values are distributed in log-scale. Since all values are in [0, 1], the log-scale is more appropriate to discriminate the values of the cost function across the state space.

Nevertheless, since in most simulations the best code is found “by chance” during the
<figure><img src="image_18.png" /><figcaption>Figure 13: Performance comparison between the best codes found in SA with the plain acceptance probability function against the initial codes in each family.</figcaption></figure>

<table><caption>Table 3: Parameters for the Simulated Annealing on each code family. The best $\beta$ parameters among the two sets of experiments are highlighted in bold.</caption><thead><tr><th>Code family</th><th>[[625, 25]]</th><th>[[1225, 65]]</th><th>[[1600, 64]]</th><th>[[2025, 81]]</th></tr></thead><tbody><tr><td>t<sub>max</sub></td><td>2400</td><td>900</td><td>450</td><td>180</td></tr><tr><td>Best β | plain-acceptance</td><td>200</td><td>10</td><td>50</td><td>200</td></tr><tr><td>Best β | log-acceptance</td><td><b>4</b></td><td><b>10</b></td><td><b>10</b></td><td><b>1</b></td></tr><tr><td>Best performance gain</td><td>5.18</td><td>273.00</td><td>77.20</td><td>2.83</td></tr></tbody></table>

intermediate iterations of the algorithm, the viability of the SA algorithm is still strongly due to the exploratory behavior it enables. Still, comparing tables 2 and 3, with roughly the same budget in terms of number of states to evaluate, Simulated Annealing managed to produce much better performance gains in comparison to plain exploration.

## 3.3 Projective Simulation

For all the experiments with Projective Simulation, we set a maximum budget for each code family in terms of the number of episodes in the training, $N$, and the maximum number of steps per episode, $L$. We also adopted a termination criterion for the episodes based on the first emission of a nonzero reward, similar to what is done by [23]. We recall that the reward function in PS is an encapsulation of the cost function, and just as [23], we adopted a sparse reward function based on a threshold $\theta$, i.e., the reward is $\lambda = 1$ if $\hat{p}_L < \theta$, and $\lambda = 0$ otherwise.

We carried out two sets of experiments, by varying the threshold for the reward function. For each code family, two values for $\theta$ were adopted, $\theta_{\text{easy}}$ and $\theta_{\text{hard}}$. Inside each of these two experiment sets, we varied the parameters $\beta$ (for the softmax function), $\gamma$ (forgetting mechanism) and $\eta$ (glow damping) over 3 values each, thus producing 27 variations. Among these 27 variations for each threshold level and for each code family, we selected only the best combination, in terms of the performance of the best code found during the training. The relation of the parameters $N$, $L$, $\theta_{\text{hard}}$, $\theta_{\text{easy}}$, the best values for $\beta$, $\gamma$ and $\eta$, as well as the best performance gains are collected in table 4.

The candidate values for the parameters $\beta$, $\eta$ and $\gamma$ were chosen in the following manner.
## 3. RESULTS AND DISCUSSION

<figure><img src="image_19.png" /><figcaption>(a) Evolution of the code performances during SA in the code family [[625, 25]].</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(b) Evolution of the code performances during SA in the code family [[1225, 65]].</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(c) Evolution of the code performances during SA in the code family [[1600, 64]].</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(d) Evolution of the code performances during SA in the code family [[2025, 81]].</figcaption></figure>

Figure 14: Evolution of the code performances throughout the execution of Simulated Annealing in each code family, for different values of $\beta$, with an acceptance probability function in terms of log-variation of the cost function.

First, the $\beta$ value is selected depending on the size of the action space, so as to yield a conversion between $h$ weights and a probability distribution where a single nonzero entry is converted to a probability of $2/3$, $1/2$ or $1/3$. Next, the $\eta$ parameter is chosen in function of $\beta$, the size of the action space and the maximum length of an episode. Three values are chosen in terms of what would be the probability that an action that received a reward once will be picked again after an entire episode. The same heuristic is followed to define the candidate values of $\gamma$.

By comparing the tables 2, 3 and 4, we see that Projective Simulation yields results only slightly better than the plain exploration, with a similar budget for each training simulation, but pales in comparison to Simulated Annealing.

In figures 16 and 17 we track the cost function of the states visited during the PS training with easy thresholds and hard thresholds, respectively. By analyzing figure 16, we see that for the three first code families, the simulation reaches a point where the episodes are terminated in just one step. This indicates that the agent has essentially learned how to achieve the desired
<figure><img src="image_23.png" /><figcaption>Figure 15: Performance comparison between the best codes found in SA with the log-scale acceptance probability function against the initial codes in each family.</figcaption></figure>

<table><caption>Table 4: Parameters for the Projective Simulation training on each code family, for the two experiment sets (hard threshold and threshold).</caption><thead><tr><th>Code family</th><th>[[625, 25]]</th><th>[[1225, 65]]</th><th>[[1600, 64]]</th><th>[[2025, 81]]</th></tr></thead><tbody><tr><td>Number of episodes</td><td>20</td><td>10</td><td>8</td><td>5</td></tr><tr><td>Max. Steps per episode</td><td>120</td><td>70</td><td>50</td><td>35</td></tr><tr><td colspan="5" style="text-align: center;">Hard threshold</td></tr><tr><td>θ<sub>hard</sub></td><td>10<sup>-2</sup></td><td>10<sup>-3</sup></td><td>2 · 10<sup>-3</sup></td><td>3 · 10<sup>-2</sup></td></tr><tr><td>Best performance gain</td><td><b>4.68</b></td><td><b>78</b></td><td><b>47</b></td><td>2.3</td></tr><tr><td>Best β</td><td>6.79</td><td>8.85</td><td>9.12</td><td>7.97</td></tr><tr><td>Best γ</td><td>4.56 · 10<sup>-4</sup></td><td>5.61 · 10<sup>-4</sup></td><td>1.72 · 10<sup>-3</sup></td><td>2.86 · 10<sup>-3</sup></td></tr><tr><td>Best η</td><td>1.90 · 10<sup>-3</sup></td><td>2.43 · 10<sup>-3</sup></td><td>1.58 · 10<sup>-3</sup></td><td>2.60 · 10<sup>-3</sup></td></tr><tr><td colspan="5" style="text-align: center;">Easy threshold</td></tr><tr><td>θ<sub>easy</sub></td><td>2 · 10<sup>-2</sup></td><td>3 · 10<sup>-3</sup></td><td>4 · 10<sup>-3</sup></td><td>4 · 10<sup>-2</sup></td></tr><tr><td>Best performance gain</td><td>4.26</td><td>54</td><td>41.9</td><td><b>2.53</b></td></tr><tr><td>Best β</td><td>6.79</td><td>8.16</td><td>8.43</td><td>8.66</td></tr><tr><td>Best γ</td><td>8.98 · 10<sup>-4</sup></td><td>0</td><td>0</td><td>0</td></tr><tr><td>Best η</td><td>0</td><td>0</td><td>0</td><td>0</td></tr></tbody></table>

threshold with ease, without exploring much the state space. This was somewhat expected for the case of easy thresholds, i.e., the best performance found becomes limited by the threshold. In figure 17, the early termination of the episodes is much less frequent. In fact, only a handful of episodes are terminated earlier than the limit, indicating that the thresholds proposed are indeed hard to attain. The limitation arising from this kind of situation is that, with few occurrences of an early termination, the agent doesn't get enough rewards, so there is not enough feedback to learn from.

## 3.4 Critique of the methods, implementations and results

Throughout the work, a number of technical limitations were presented and briefly discussed. Some of these limitations hint at potential directions to be explored in further work.
3. RESULTS AND DISCUSSION

<figure><img src="image_24.png" /><figcaption>Figure 16: Evolution of the cost function of the states visited throughout the PS training, for a reward function set with an easy threshold.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Figure 17: Evolution of the cost function of the states visited throughout the PS training, for a reward function set with a hard threshold.</figcaption></figure>
* **Importance Sampling:** It was pointed out in section 2.1.3 that the plain Monte Carlo estimation had a limitation in its precision related to the difficulty to estimate the probability of rare events. In the literature of Monte Carlo methods, importance sampling is one of the techniques that may help overcome this problem. We tried to apply such techniques in the case of estimating failure rates for very low noise levels, without success. The difficulty related to this problem is that in order to produce good quality estimations via importance sampling, the sampling distribution must be as close as possible to that of the random variable whose expected variable is being estimated. In the case of quantum codes, it is not straightforward to understand how to tune this distribution. It is known that in the low noise regime, the failure rate is dominated by logical errors of very small weight, and thus it is not sufficient to simply increase the erasure rate to make the failures more frequent, because that will lead to failures caused by errors of large weight, causing the failure rate to be underestimated.

* **Graph Canonization:** In section 2.1.1, we briefly mentioned that many different states may describe the same code, most notably those states whose graphs are isomorphic. When discussing the exploration algorithm, we mentioned the work of [42], that describes how the number of states in our space scales with the numbers of bits and checks. While it is not clear how much the space gets reduced when one quotients the space by some equivalence relation, such as grouping together isomorphic graphs, it is nevertheless a way to reduce the overall number of states, and also to allow for the caching of the cost function for equivalent states, potentially speeding up the simulations. However, my attempts at implementing the identification of equivalent states stumbled upon impractical computational costs associated to the graph isomorphism problem, for which no polynomial algorithm is known.

* **Decoder complexity and exploration budget:** Due to the $O(N^3)$ complexity of the performance estimation procedure, even a C/C++ implementation of the cost functions falls short of efficiency, becoming thus the main computational bottleneck of our approach. In effect, this elevated costs hinders the applicability of exploration methods such as Projective Simulation, which, like many RL algorithms, tends to be very sample inefficient. A possibility for further investigation would be to train a Graph Neural Network (GNN) to learn the cost function from annotated data and replace the Monte Carlo simulation by the GNN's predictions.

* **Parameter tuning and reward engineering in PS:** RL algorithms are known for being hard to tune, and the strategies of reward engineering are applicable on a case-by-case basis. In the current work, I believe that alternative reward functions or episode termination criteria could be more effective than the one adopted here. However, due to time and computational budget constraints, not many experiments were carried out with PS. As a possibility for further investigation, we can replace PS with other, more standard RL algorithms, particularly actor-critic algorithms. By employing a pair of GNNs to do the action proposal and the reward estimation, we could potentially reduce the cost of evaluating a state and also overcome one of the limitations of projective simulation, which is the lack of generalization of its agent's experience to unvisited states.
Conclusion

In this work, we presented a brief background on Quantum Error Correction, and addressed the problem of optimizing quantum LDPC codes, in particular through the construction of hypergraph product codes. These codes are interesting because they provide a large source of qLDPC codes, and also for allowing the formulation of the optimization problem in terms of the classical codes making up the product.

In section 2 we presented our formulation of the optimization problem in terms of a search problem in a space of codes, that can be traversed by the application of local deformations of the codes. The metric of performance chosen was the estimation of the logical error rate of the quantum code under the quantum erasure channel. We introduced the Projective Simulation framework, as well an alternative optimization method, the Simulated Annealing algorithm.

We carried out experiments testing each of the proposed optimization methods, effectively obtaining better codes for four families of qLDPC codes of interest. An aspect worth pointing out is that codes that were optimized in a particular noise level also displayed better performances in broader ranges of noise levels. While both methods were capable to produce better codes than other standard, ad hoc techniques, such as the edge progressive growth algorithm, the results showed that the effect of random exploration of the state spaces is the main responsible for the gains in performance. Nevertheless, the tuning of the Simulated Annealing algorithm was shown to be much more successful than that of the Projective Simulation approach.

The analysis of the particular structure of the codes achieving better performances can be an interesting direction of investigation for future work, both from the mathematical, analytical point of view and as a potential source of heuristic methods for optimizing these codes.
References

[1] Ville Lahtinen and Jiannis Pachos. A short introduction to topological quantum computation. *SciPost Physics*, 3(3), September 2017.

[2] Diego Ruiz, Jérémie Guillaud, Anthony Leverrier, Mazyar Mirrahimi, and Christophe Vuillot. Ldpc-cat codes for low-overhead quantum computing in 2d, 2024.

[3] Kosuke Fukui. High-threshold fault-tolerant quantum computation with the gottesman-kitaev-preskill qubit under noise in an optical setup. *Physical Review A*, 107(5), May 2023.

[4] Nicholas Connolly, Vivien Londe, Anthony Leverrier, and Nicolas Delfosse. Fast erasure decoder for hypergraph product codes. *Quantum*, 8:1450, August 2024.

[5] Nicolas Delfosse and Naomi H. Nickerson. Almost-linear time decoding algorithm for topological codes. *Quantum*, 5:595, December 2021.

[6] Nicolas Delfosse and Gilles Zémor. Linear-time maximum likelihood decoding of surface codes over the quantum erasure channel. *Phys. Rev. Res.*, 2:033042, Jul 2020.

[7] Joschka Roffe, David R. White, Simon Burton, and Earl Campbell. Decoding across the quantum low-density parity-check code landscape. *Physical Review Research*, 2(4), December 2020.

[8] Laia Domingo Colomer, Michalis Skotiniotis, and Ramon Muñoz-Tapia. Reinforcement learning for optimal error correction of toric codes. *Physics Letters A*, 384(17):126353, June 2020.

[9] Ye-Hua Liu and David Poulin. Neural belief-propagation decoders for quantum error-correcting codes. *Physical Review Letters*, 122(20), May 2019.

[10] Philip Andreasson, Joel Johansson, Simon Liljestrand, and Mats Granath. Quantum error correction for the toric code using deep reinforcement learning. *Quantum*, 3:183, September 2019.

[11] Nikolas P. Breuckmann and Jens Niklas Eberhardt. Quantum low-density parity-check codes. *PRX Quantum*, 2(4), October 2021.

[12] Omar Fawzi, Antoine Grospellier, and Anthony Leverrier. Constant overhead quantum fault-tolerance with quantum expander codes. In *2018 IEEE 59th Annual Symposium on Foundations of Computer Science (FOCS)*. IEEE, October 2018.

[13] Jean-Pierre Tillich and Gilles Zemor. Quantum ldpc codes with positive rate and minimum distance proportional to the square root of the blocklength. *IEEE Transactions on Information Theory*, 60(2):1193–1202, February 2014.
[14] Dolev Bluvstein, Simon J. Evered, Alexandra A. Geim, Sophie H. Li, Hengyun Zhou, Tom Manovitz, Sepehr Ebadi, Madelyn Cain, Marcin Kalinowski, Dominik Hangleiter, J. Pablo Bonilla Ataides, Nishad Maskara, Iris Cong, Xun Gao, Pedro Sales Rodriguez, Thomas Karolyshyn, Giulia Semeghini, Michael J. Gullans, Markus Greiner, Vladan Vuletić, and Mikhail D. Lukin. Logical quantum processor based on reconfigurable atom arrays. *Nature*, 626(7997):58–65, December 2023.

[15] Pavel Panteleev and Gleb Kalachev. Asymptotically good quantum and locally testable classical ldpc codes, 2022.

[16] Anthony Leverrier and Gilles Zémor. Quantum tanner codes, 2022.

[17] Xiao-Yu Hu, Evangelos Eleftheriou, and Dieter-Michael Arnold. Progressive edge-growth tanner graphs. *GLOBECOM'01. IEEE Global Telecommunications Conference (Cat. No.01CH37270)*, 2:995–1001 vol.2, 2001.

[18] Jan Olle, Remmy Zen, Matteo Puviani, and Florian Marquardt. Simultaneous discovery of quantum error correction codes and encoders with a noise-aware reinforcement learning agent, 2024.

[19] Vincent Paul Su, ChunJun Cao, Hong-Ye Hu, Yariv Yanay, Charles Tahan, and Brian Swingle. Discovery of optimal quantum error correcting codes via reinforcement learning, 2023.

[20] Caroline Mauron, Terry Farrelly, and Thomas M. Stace. Optimization of tensor network codes with reinforcement learning, 2023.

[21] Hans J. Briegel and Gemma De las Cuevas. Projective simulation for artificial intelligence, 2013.

[22] Alexey A. Melnikov, Adi Makmal, Vedran Dunjko, and Hans J. Briegel. Projective simulation with generalization. *Scientific Reports*, 7(1), October 2017.

[23] Hendrik Poulsen Nautrup, Nicolas Delfosse, Vedran Dunjko, Hans J. Briegel, and Nicolai Friis. Optimizing Quantum Error Correction Codes with Reinforcement Learning. *Quantum*, 3:215, December 2019.

[24] Austin G. Fowler, Matteo Mariantoni, John M. Martinis, and Andrew N. Cleland. Surface codes: Towards practical large-scale quantum computation. *Physical Review A*, 86(3), September 2012.

[25] Nicolas Delfosse, Pavithran Iyer, and David Poulin. A linear-time benchmarking tool for generalized surface codes, 2016.

[26] Antoine Grospellier. *Constant time decoding of quantum expander codes and application to fault-tolerant quantum computation*. Theses, Sorbonne Université, November 2019.
[27] M.A. Nielsen and I.L. Chuang. *Quantum Computation and Quantum Information: 10th Anniversary Edition*. Cambridge University Press, 2010.

[28] R. Gallager. Low-density parity-check codes. *IRE Transactions on Information Theory*, 8(1):21–28, 1962.

[29] Claude E. Shannon. A mathematical theory of communication. *Bell System Technical Journal*, 27:379–423, 623–656, 1948.

[30] Daniel Gottesman. Stabilizer codes and quantum error correction. *Ph. D. Thesis*, 1997.

[31] A. R. Calderbank and Peter W. Shor. Good quantum error-correcting codes exist. *Phys. Rev. A*, 54:1098–1105, Aug 1996.

[32] Andrew Steane. Multiple-particle interference and quantum error correction. *Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences*, 452(1954):2551–2577, 1996.

[33] Yue Wu, Shimon Kolkowitz, Shruti Puri, and Jeff D. Thompson. Erasure conversion for fault-tolerant quantum computing in alkaline earth rydberg atom arrays. *Nature Communications*, 13(1), August 2022.

[34] Nicolas Delfosse and Gilles Zémor. Upper bounds on the rate of low density stabilizer codes for the quantum erasure channel, 2012.

[35] Victor V. Albert and Philippe Faist. Quantum low-density parity-check (qldpc) code. In *The Error Correction Zoo*. https://errorcorrectionzoo.org/c/qldpc, 2022.

[36] Dimiter Ostrev, Davide Orsucci, Francisco Lázaro, and Balazs Matuz. Classical product code constructions for quantum Calderbank-Shor-Steane codes. *Quantum*, 8:1420, July 2024.

[37] Martin Albrecht and Gregory Bard. *The M4RI Library*. The M4RI Team. Available at https://bitbucket.org/malb/m4ri.

[38] Alexey A. Melnikov. *Projective Simulation*. Available at https://github.com/qic-ibk/projectivesimulation.

[39] Aric Hagberg, Pieter Swart, and Daniel S Chult. Exploring network structure, dynamics, and function using networkx. Technical report, Los Alamos National Lab.(LANL), Los Alamos, NM (United States), 2008.

[40] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Haldane, Jaime Fernández del Río, Mark Wiebe, Pearu Peterson, Pierre Gérard-Marchant, Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi, Christoph Gohlke, and Travis E. Oliphant. Array programming with NumPy. *Nature*, 585(7825):357–362, September 2020.
[41] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Sté-fan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay May-orov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimr-man, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. *Nature Methods*, 17:261–272, 2020.

[42] Jeffrey W. Miller and Matthew T. Harrison. Exact sampling and counting for fixed-margin matrices. *The Annals of Statistics*, 41(3), June 2013.

[43] Dimitris Bertsimas and John Tsitsiklis. Simulated annealing. *Statistical Science*, 8(1):10–15, 1993.