Degenerate Quantum LDPC Codes With Good Finite Length Performance

We study the performance of medium-length quantum LDPC (QLDPC) codes in the depolarizing channel. Only degenerate codes with the maximal stabilizer weight much smaller than their minimum distance are considered. It is shown that with the help of OSD-like post-processing the performance of the standard belief propagation (BP) decoder on many QLDPC codes can be improved by several orders of magnitude. Using this new BP-OSD decoder we study the performance of several known classes of degenerate QLDPC codes including hypergraph product codes, hyperbicycle codes, homological product codes, and Haah's cubic codes. We also construct several interesting examples of short generalized bicycle codes. Some of them have an additional property that their syndromes are protected by small BCH codes, which may be useful for the fault-tolerant syndrome measurement. We also propose a new large family of QLDPC codes that contains the class of hypergraph product codes, where one of the used parity-check matrices is square. It is shown that in some cases such codes have better performance than hypergraph product codes. Finally, we demonstrate that the performance of the proposed BP-OSD decoder for some of the constructed codes is better than for a relatively large surface code decoded by a near-optimal decoder.


Introduction
Quantum error-correcting codes are considered as an essential component in the current architectures of quantum computers due to the inherently faulty nature of the quantum hardware. Topological quantum codes [1,2,3] are among the quantum codes with the highest known noise thresholds. These codes have sparse parity-check matrices and thus belong to the class of quantum LDPC (QLDPC) codes. Moreover, they are highly degenerate, which means that the minimum distance is much higher than the weight of the stabilizers. It is important that these codes also have decoding algorithms with near to optimal performance [1,4,5]. Though the thresholds of topological codes are relatively high, their dimensions are usually much smaller than for general QLDPC codes of the same length (it is constant for the surface and color codes). There are also interesting classes of topological quantum codes with very large dimensions (e.g., hyperbolic codes [6] have a constant rate). However, such codes usually have much smaller minimum distances compared to the surface and color codes.
Recently there have been proposed a number of interesting families of degenerate QLDPC codes (e.g., hypergraph product codes [7] and homological product codes [8]) with very good asymptotic parameters. Nevertheless their practical error-correcting performance for relatively small code lengths (n < 1000) is largely unexplored, and it is not clear whether their performance is competitive to the best known topological codes. From our point of view, the difficulty of constructing degenerate QLDPC codes with good practical performance is mostly related to the following two issues.
1. Asymptotically good constructions may not necessarily produce the best QLDPC codes for relatively small code lengths. Indeed, in [9,10] the construction of hypergraph product codes was further improved and generalized. Although the asymptotic characteristics of the improved codes are the same as before, their parameters such as the rate and the minimum distance are much better for smaller lengths.
2. The performance of the known QLDPC codes is far from optimal under the state-of-the-art decoders (including the binary and non-binary BP decoders, and their modifications). The performance degradation is usually attributed [11] to the unavoidable 4-cycles in the corresponding Tanner graphs and to a large number of degenerate errors. While the number of 4-cycles can be significantly reduced by using CSS codes [12] without 4-cycles in the paritycheck matrices H X and H Z , the number of degenerate errors can not be easily reduced for highlydegenerate codes.
In this paper, we try to address both mentioned issues. In the first part of the paper, we introduce a new enhancement of the standard BP decoder for QLDPC codes (both binary and non-binary versions are allowed) using a variant of the well-known decoding algorithm for short classical codes called the ordered statistics decoding (OSD) [13,14]. This new post-processing algorithm works only if the BP decoder fails to find a recovery Pauli operator that gives the correct syndrome. We suppose here that the QLDPC parity-check matrix H is represented in the binary form 1 . The algorithm starts from finding a reliable information set [15] for H based on the soft decisions obtained by the BP decoder. Then it makes hard decisions for the bits from this information set and flips the w most unreliable of them in order to find the 2 w corresponding recovery Pauli operators that return the corrupted quantum state to the coding space. Finally, it selects a recovery operator with the minimum weight and applies it to the corrupted quantum state. Thus this new combined BP-OSD decoder, in contrast to the standard BP, always returns the recovery operator that moves the corrupted codeword back to the code space. We show that with the help of this OSD-like post-processing the performance of the standard BP decoder on many degenerate QLDPC codes can be improved by several orders of magnitude. As we can see from its brief description above, it is not restricted to CSS codes and can also be used in conjunction with any decoder, different from BP, that provides soft decisions.
In the second part of the paper, we construct a number of new relatively small codes using two families of degenerate codes: the generalized bicycle (GB ) codes introduced in [10] and a new large family 2 of QLDPC codes introduced in this paper. We call this new family generalized hypergraph product (GHP ) codes. It contains the GB codes and the class of hypergraph product (HP) codes [7], where one of the two parity-check matrices used in the product is square. We also derive new formulas for the dimension of GB and GHP codes and show how to design codes of high dimension with good error correction performance. It is interesting to note that some of the constructed GB codes have an additional property that their syndromes are protected by small BCH codes, which may be useful for the faulttolerant syndrome measurement.
We also study the performance of the proposed BP- 1 For a QLDPC code on n qubits with m stabilizers the matrix H is an m × 2n binary matrix. 2 Since the first version of the current work was released, the codes from this family have been shown to have very large minimal distances [16,17,18] and found interesting applications in geometry [19]. In [17] they were further generalized and called lifted product codes.
OSD decoder on many known classes of degenerate QLDPC codes, including the already mentioned hypergraph product codes, hyperbicycle codes [9,10], the homological product codes [8], and Haah's cubic codes [20]. We compare their performance with the performance of the codes, constructed in this work. We show that in many cases the new codes with similar performance have better parameters such as the code length and the rate. Besides that, we compare the new BP-OSD decoder with other known modifications of the BP decoder such as the random perturbation [21], the enhanced feedback [22], and the matrix augmentation [23] algorithms. We compare all the above mentioned algorithms on the new [[1270, 28]] code from the class of GHP codes mentioned earlier 3 . We show that the performance of the new decoding algorithm on this code is significantly better. Moreover, we also show that the performance of this code under the BP-OSD is even better than the performance of the [[1201, 1,25]] surface code under the near-optimal MPS-based decoder proposed in [5]. We also demonstrate that the performance under BP-OSD of one of Haah's cubic codes, which have local stabilizers in 3D, is also very good.
The remainder of the paper is organized as follows. Section 2 contains some background material, where we review standard definitions and fix notations. In Section 3, we describe the BP-OSD algorithm and compare its performance with other known modifications of BP. In Section 4, we study GB codes and construct some new codes with good performance. In Section 5, we introduce and study a new large family of GHP codes and show that it generalizes the class of hypergraph product codes in the case when one of the parity-check matrices is square. In the last section, we give some final remarks. The paper also contains three appendixes. Appendix A has some supplementary material on the ring of circulants. Appendix B contains a description of all the codes used in the simulations. Finally, we show some additional simulations in Appendix C.

Basic facts and definitions
In this section, we fix notations and briefly recall some standard definitions related to classical and quantum LDPC codes. See [11] for a good review of these topics.

Classical codes
Let n be a natural number. In what follows, we denote by [n] the set {1, . . . , n}. Consider a finite field F q 3 Since this GHP code is quasi-cyclic, it can also be obtained (up to some qubit permutations) as a special case of the hyperbicycle code construction. and an n-dimensional vector space F n q over F q . A linear [n, k] q code is a k-dimensional subspace C ⊆ F n q , where the parameters n and k are called the length and the dimension of C, respectively. We denote the dimension k of the code C by dim C. The rate of the code C is equal to k/n. The elements of C are called codewords. The Hamming distance d(v, v ) between vectors v, v ∈ F n q is the number of positions in which they differ. The parameter is called the minimal distance of C. By definition, we put d(C) = ∞ when k = 0. It is easy to see that d(C) is equal to the minimal weight |c| of non-zero codewords, where the weight |c| is the number of non-zero components in c. When d(C) = d for a linear [n, k] q code C, we say that C is an [n, k, d] q code. A linear [n, k, d] q code is usually defined either as the row space of a matrix G called the generator matrix or as the kernel of a matrix H called the parity-check matrix. It is easy to see that GH T = 0, rk G = k, and rk H = n − k. In what follows we mostly consider binary linear codes (i.e., q = 2), in which case we use a shorter notation: [n, k] or [n, k, d].

Quantum stabilizer codes
The quantum analogs of classical linear codes are quantum stabilizer codes introduced in [24]. To define them we need a number of supporting definitions. Consider an n-qubit Hilbert space C 2 n = (C 2 ) ⊗n . A Pauli operator on n qubits is an operator P = αP 1 ⊗· · ·⊗P n on the space C 2 n , where α ∈ {±1, ±i}, P i ∈ {I, X, Y, Z}. Here I is the identity and X, Y, Z are the Pauli 2 × 2 matrices. The weight wt (P ) of a Pauli operator P is the number of non-identity components in the tensor product. The set P n of all Pauli operators P on n qubits is a non-commutative group under the operator multiplication called the n-qubit Pauli group. While the group P n is non-commutative, if we forget about the global phase factors α of Pauli operators we obtain the quotient group P n /{±I, ±iI}, where I is the identity operator on C 2 n . This quotient group is isomorphic to the commutative group Z 2n 2 , and the isomorphism is given by the following map: (1) It is well known that two n-qubit Pauli operators P and P commute if f for their binary representations (x|z), (x |z ) ∈ F 2n 2 we have where a, b = i a i b i is the dot product in F n 2 .
Denote by P * n the subset of Pauli operators P ∈ P n with the global phase factor α = 1. A stabilizer group S is a commutative subgroup of the Pauli group P n such that −I ∈ S. The group S is usually generated by m Pauli operators S 1 , . . . , S m ∈ P * n called stabilizers generators, i.e., S = S 1 , . . . , S m . We say that S 1 , . . . , S m are independent if none of them can be obtained (up to a global phase factor α) from the others by the group multiplication.
Let us recall that a quantum stabilizer [[n, k, d]] code is a 2 k -dimensional subspace of the n-qubit space C 2 n defined as the common +1-eigenspace for a set of m stabilizers S 1 , . . . , S m ∈ P * n : where the group S(C) = S 1 , . . . , S m is called the stabilizer group of C. It can be shown that m n − k. But if the stabilizers are independent, we get m = n − k.
Here the parameter d = d(C) is called the minimal distance 4 of C and is equal to the minimal possible weight of a Pauli operator P ∈ P * n that commutes with all the stabilizers S 1 , . . . , S m but P ∈ S(C).
A Pauli operator P ∈ P * n is usually interpreted as an error operator (called a Pauli error ) that can corrupt a quantum system and cause it to go from a state |ψ to P |ψ . However, for every quantum stabilizer code C, it is not hard to show that P |ψ = |ψ for all |ψ ∈ C if f P ∈ S(C). Hence we see that in this case, not all Pauli errors P ∈ P * n can harm the state |ψ . We call a Pauli error P ∈ P * n degenerate for a code C if P |ψ = |ψ for all |ψ ∈ C and non-degenerate otherwise. We see that the elements from P ∈ S(C) are precisely the degenerate Pauli errors and the minimum distance d(C) is the minimal possible weight of a non-degenerate Pauli error. A stabilizer code C is called degenerate if it has degenerate Pauli errors P of weight wt (P ) < d(C).
If we apply the binary mapping (1) to the stabilizer generators S 1 , . . . , S m of an [[n, k, d]] code C, we obtain the m × 2n binary matrix called the parity-check matrix of C, where each row corresponds to a stabilizer generator. We do not require for the matrix H to be full rank (i.e, rk H = n−k m). Using (2) it is not hard to see that the null space of H is exactly the set of all vectors (z, x) ∈ F 2n 2 such that (x, z) is the binary mapping of a Pauli error P ∈ P * n that commutes with all the stabilizer generators S 1 , . . . , S m .
We see that the matrix H = (H X | H Z ) is not an arbitrary binary m × 2n matrix since the stabilizer generators S 1 , . . . , S m corresponding to its rows should com-mute with each other. From equation (2) it easily follows that this restriction can be formulated as the following commutativity condition: A very important subclass of quantum stabilizer codes is Calderbank-Shor-Steane (CSS) codes introduced in [25,26]. A quantum CSS code is a stabilizer code, where non-identity components in the tensor product of each stabilizer generator are either all equal to X or all equal to Z. Hence the parity-check matrix H of a CSS code of length n can be represented in the following form: where H X , H Z are binary matrices with n columns. In this special case, commutativity condition (4) can be rewritten as: We can easily verify that the dimension k of the corresponding CSS code of length n is given by the formula:

Classical and quantum LDPC codes
A classical low density parity check (LDPC) code [27] is a linear code defined by a sparse binary parity-check matrix H = (h ij ) m×n . The sparseness usually means that the weights of all rows and columns in H are upper bounded by some universal constant as the code length n grows in an infinite family of codes. When we consider an LDPC code defined by a paritycheck matrix H, it is helpful to define the bipartite graph T = T (H) called the Tanner graph [28]. In this graph the first part of nodes v 1 , . . . , v n (called the vnodes) corresponds to the columns of H (the variables), the second part of nodes c 1 , . . . , c m (called the c-nodes) corresponds to the rows of H (the checks), and we con- . If the parity-check matrix H is (w c , w r )-regular (i.e., each column has weight w c and each row has weight w r ), then the corresponding Tanner graph is also (w c , w r )-regular (i.e., each v-node has degree w c and each c-node has degree w r ). We say that an LDPC code is w-limited if the degree of each node in its Tanner graph is upper bounded by w. It is obvious that any LDPC code with (w c , w r )-regular parity-check matrix is max(w c , w r )-limited.
There are a number of decoding algorithms for classical LDPC codes, but the most frequently used one is the belief propagation (BP) decoder [27], also known as the message passing decoder or the sum-product decoder [29]. It assigns the a priori probability distributions of individual bits in the codeword (obtained from the channel) to the v-nodes of the Tanner graph and iteratively updates the posterior probability distributions for each bit. Once some maximal iteration number limit is reached, the decoder combines the a priori and the calculated posterior probability distributions (the soft decisions) to produce the optimal binary decision (called the hard decision) for each individual bit.
An important property of a parity-check matrix H defining an LDPC code is the girth of the corresponding Tanner graph T = T (H), which is equal to the length of the shortest cycle in T . It is well known that short cycles in the Tanner graph degrade the performance of the BP decoder. At the same time, it was observed that LDPC codes without 4-cycles (i.e., when the girth is at least 6) perform very well in practice. However, this practical observation is not fully investigated from theoretical point of view.
A quantum LDPC (QLDPC) is a stabilizer [[n, k, d]] code with a sparse parity-check matrix H. We can also introduce the Tanner graph T = T (S) for any stabilizer [[n, k, d]] code C defined by a set of stabilizer generators S = {S 1 , . . . , S m }. In the case of stabilizer codes, the v-nodes correspond to n qubits and the c-nodes to the stabilizer generators S 1 , . . . , S m , and we connect a cnode with a v-node if the corresponding stabilizer acts nontrivially on the corresponding qubit. As in the case of classical LDPC codes, we say that a QLDPC code is w-limited if the degree of each node in its Tanner graph is upper bounded by w. This property is much more important in the quantum case due to the faulty nature of the current quantum hardware. It is clear that any CSS code with (w c , w r )-regular matrices H X and H Z is max(2w c , w r )-limited.

OSD-like post-processing for BP
In this section, we describe a new OSD-like postprocessing algorithm that can be used after the BP decoder for QLDPC codes. Before we give its detailed description we consider two simple modifications of the OSD decoder for classical linear codes. These modifications will be used as the main components in the OSD-like post-processing algorithm for quantum codes. We should warn the reader that these modifications of the standard OSD decoder are not intended to improve its performance for classic LDPC codes. We introduce them because these algorithms eventually will be used in the OSD post-processing algorithm (called qOSD) for quantum codes described in Section 3.2. For example, one of the main differences between classical and quantum codes is that we have to use the syndrome decoder in the quantum case. Hence we consider only syndrome OSD decoders here since we are going to use them as components of the decoder for quantum codes.

Syndrome OSD post-processing algorithm
Starting from this section we will often use the following notations: • If M is a matrix, then M i denotes its i-th column.
and π ∈ S n is a permutation, then for every vector v = (v 1 , . . . , v n ) and matrix M = (M 1 , . . . , M n ) we define: is an index set, then the index setĪ = [n] \ I is called its complement.

Now let us recall the definition of an information set [30] of a code, which has an important role in the OSD decoding. A set of indices
Clearly, I is an information set of C if f |I| = k, and for any two codewords c, c ∈ C such that c I = c I we always have c = c . Therefore for any information set I the corresponding k bits can be used to uniquely recover any vector v ∈ F n 2 if we also know its syndrome s = Hv. Hence we can consider the corresponding encoding map E s I : F k 2 → F n 2 such that for any u, v, and s we have: . It is easy to verify that E 0 I is the systematic encoding map for the code C where the information bits u are at the positions corresponding to I. Remark. If G is a generator matrix, and H is a paritycheck matrix of a linear code C; then it can be shown denotes the family of indices I such that the collection of columns {M i } i∈I is a basis for the column space of M ; then the family of all information sets of C coincides with B(G), while the family of their complements coincides with B(H). Since the set of all linearly independent columns of any matrix gives us a matroid 5 , 5 A matroid is defined by a non-empty collection I of subsets (called the independent sets) from some set E (called the ground set) such that: (1) I is closed under taking subsets; (2) for any two A, B ∈ I if |A| < |B| then A ∪ {b} ∈ I for some b ∈ B \ A. Any maximal (by inclusion) independent set is called a basis. then for any positive real numbers w 1 , . . . , w k we can efficiently find arg max in a greedy fashion [31,32]. Here in the right part we The main idea of the syndrome OSD decoder is as follows. Consider a linear [n, k] code C defined by a paritycheck matrix H. Let c = c + e be a corrupted version of a codeword c ∈ C, where e ∈ F n 2 is the corresponding random error vector. Given the syndrome s = Hc = He we want to find the error vector e and thus recover the codeword c = c − e. Suppose that in addition we are given an estimate 6 of the error probability p i = P(e i = 1) for each i ∈ [n]. From these estimates our best guess of the error vector e would be the hard decisions vectorê, whereê i = 1 when p i > 1/2, andê i = 0 otherwise 7 . However, when we only use the error probability estimates, we often have that Hê = s, and thusê = e. Nevertheless, even in such cases, some components ofê are equal to e, and we can still try to useê to find e if we also take into account that Therefore, to recover e, one may try to traverse through different information sets I in the hope that for one of them we haveê I = e I , and thus e = E s I (ê I ).
Unfortunately we do not know for which indices i ∈ [n] we haveê i = e i . Hence it makes sense to find an information set I with the indices that are as reliable as possible, where the reliability of an index i ∈ [n] is the probability If we assume that the components in the random error vector e are mutually independent, then it follows that the probability of successful decoding P(e = E s I (ê I )) for I is equal to ρ(I) = i∈I ρ i . Thus if we want to maximize this probability we need the most reliable information set I, i.e., the one 8 with the maximal possible value of ρ(I). Finding such a set may at first sight look like a prohibitively hard task. However, if instead of ρ(I) we consider ln ρ(I) = i∈I w i , where w i = ln ρ i ; then Algorithm 1: Syndrome OSD-w algorithm Input: target weight function wt (·), binary parity-check matrix H, syndrome vector s ∈ F m 2 , vector of hard decisionsê ∈ F n 2 ; Output: error vector e ∈ F n 2 such that He = s; we can see that the most reliable information set I is given by (7), and, as we mentioned earlier, it is possible to find I using a greedy algorithm [31,32]. This greedy algorithm is the main part of the OSD decoder. Since in our case the code is defined by a parity-check matrix H, it is more convenient to use the right part of (7), which means that we want to find the index set J corresponding to the least reliable basis {H i } i∈J for the column space of H, i.e., J ∈ B(H) with the smallest ρ(J). For simplicity, we assume that the columns of H are already rearranged such that the reliability of the corresponding positions increases: ρ 1 · · · ρ n . In this case, it is not hard to see that the collection {H i } i∈J of the first r = rk H linearly independent columns is the least reliable basis [32]. We can find J by applying Gaussian elimination to equation (8), which gives us, in time O(n 3 ), the following equation:H e =s.
Then J is the index set of the first r pivot columns inH, i.e.,H J = (δ ij ) r×r . Moreover, from (7) it follows that the index set I =J is the most reliable information set, and we can find the error vector e from (10) in a very straightforward way sinceH J is the identity matrix. The described above syndrome decoding algorithm is usually called order-0 OSD decoder. As we see, it first finds the most reliable information set I and then recovers the unique vector e ∈ F n 2 subject to the following conditions: In fact, we can further improve the error-correcting performance by using order-w OSD decoder (abbreviated as OSD-w). In this modification, after we find the most reliable information set I, we look at all error vectors e obtained from equation (10) by setting the first w least Algorithm 2: Fast syndrome OSD-0 algorithm Input: binary parity-check matrix H, syndrome vector s ∈ F m 2 , vector of hard decisionsê ∈ F n 2 ; Output: error vector e ∈ F n 2 such that He = s; reliable positions from I to all possible values x ∈ F w 2 (we can obtain them using the encoding operator E s I ). Finally, we select an error vector e of the minimal weight wt (e) among all the 2 w obtained error vectors. Here the weight function wt (·) depends on the specific channel we use. In the case of depolarizing noise, this weight function is just the weight of the corresponding Pauli operator.
A simplified pseudocode of the above OSD-w algorithm is shown in Algorithm 1. Here in the two last lines, we used a shorthand notation R x I v for the result of the replacement of the subvector v I in the vector v by x ∈ F |I| 2 . Therefore, R x [w]ê I is obtained fromê I if we replace its first w bits by x ∈ F w 2 . Since we assume that the columns of H are already sorted according to the reliabilities, these first w positions in I are the least reliable ones. We should also note that the main cycle of this algorithm, which finds J (lines 1-7), can be efficiently implemented using Gaussian elimination, as already discussed above.
Remark. There are some differences between the standard order-w OSD decoder from [13,14] and the proposed OSD-w post-processing algorithm. For example, we try all 2 w bit flips of the w least reliable information bits (line 8), while in the standard OSD we try all the i w k i bit flips of no more than w bits. As we will see in Section 3.3, if the OSD-w decoder is used as a post-processing algorithm after the BP decoder for QLDPC codes, then it can radically improve the error correcting performance. In fact, in many cases even OSD-0 is enough, and thus the computational cost of the decoding is O(n 3 ). Moreover, in the case of the OSD-0, there is no need to do full Gaussian elimination for H. Let I be the most reliable information set found by the OSD-0 algorithm, then we can stop extending the index set J in Algorithm 2 when rk [H J , s ] = rk H J , and put e = R x Jê , where s = s + HJêJ , and x is the unique 9 solution of H J x = s . Indeed, we get Since we also have I ⊆J, then e satisfies conditions (11) and thus is the error vector found by the OSD-0 algorithm. This observation gives us a faster version of the OSD-0 algorithm (see Algorithm 2).

Modified OSD post-processing algorithm for stabilizer codes
In this subsection, we show how to adapt the OSD decoder from the previous subsection to a scenario where it is used as a post-processing algorithm after the BP decoder for QLDPC codes. There are a number of quantum noise models. In this paper, we consider the depolarizing channel only. However many of our ideas may be used for other memoryless quantum noise models. In the depolarizing channel model with error probability p, a quantum state |ψ ∈ C 2 n is subject to a random Pauli error where all E i are i.i.d, and for all i ∈ [n] we have: In what follows, it is convenient for our goals to represent (with a small abuse of notation and terminology) the set of matrices P * 1 = {I, X, Y, Z} by the elements of the finite field F 4 , where I is represented by 0 ∈ F 4 , and the Pauli matrices X, Y, Z by the three non-zero elements from F 4 . To distinguish these finite field elements from the corresponding matrices we denote the former using a different font: i, x, y, z. Further, we represent a Pauli vector from P * n by the corresponding vector from F n 4 , which we also call a Pauli vector. Using this conventions, we represent the binary m×2n paritycheck matrix of a stabilizer code C (see equation (3)) by the corresponding m × n matrix over F 4 and call it the stabilizer matrix of C.
Since we consider the depolarizing channel, for the best performance, it is better to use the non-binary version of the syndrome BP decoder (see [21], [11, Algorithm 1]), which also takes into account the correlations between X and Z errors in qubits. To describe the OSD algorithm in this case we need some extra notations. 9 The solution is unique since rk H J = |J|.
; then we can define the binary vectors: We also need the inverse The main motivation for these, not very standard, notations is as follows. If e ∈ F n 4 is a Pauli error, then it is easy to check that If we have the Z error in the second qubit, then we get the error vector e = (i, z, i, i, i) ∈ F 5 4 , its binary representation b(e) = (0, 0, 0, 1, 0, 0, 0, 0, 0, 0) ∈ F 10 2 , and the corresponding syndrome vector s = (0, 1, 0, 1).
Since the depolarizing channel is non-binary, we need some further adjustments in the syndrome OSD decoder (Algorithms 1 and 2) to use it as a post-processor after the non-binary BP decoder. We call this modified version the order-w qOSD algorithm, which is defined via the order-w OSD algorithm from the previous section as follows: Here we use the symplectic weight 10 |·| Sp as the target weight function wt (·) for the OSD algorithm, which is 10 It is easy to see that the symplectic weight |x| Sp is equal to the weight wt (E) of the Pauli error E ∈ P * n represented in the binary form by the vector x ∈ F 2n 2 .
defined as: As we see, the input for the qOSD algorithm includes the m × n stabilizer matrix H over F 4 , the syndrome s ∈ F m 2 , and the error vectorê ∈ F 4 , which is the hard decisions vector for the result of the non-binary BP decoder (see the next subsection).
Remark. The symplectic weight |·| Sp is used above for the depolarizing channel. For the variant of the depolarizing channel, where the X and Z errors are independent, the standard binary Hamming weight |·| should be used instead. In general, the optimal choice of the target weight function wt (·) is the one, where the most probable errors are of the least possible weight.

Main decoding algorithm (BP-OSD)
In this subsection, we describe our main decoding algorithm (see Algorithm 3) called the BP-OSD or BP-OSD-w (if we want to emphasize the OSD order w). It consists of two stages: the non-binary BP decoding (lines 1-4) and the OSD post-processing (lines 6-8).
In fact, at the first stage, any decoding algorithm may be used instead of BP if it can provide soft decisions. The goal of this stage is either to find the error vector e ∈ F n 4 or to obtain the probabilities P(e i = E), i ∈ [n], for each type of the Pauli error E ∈ {i, x, y, z}.
First, we run the BP decoder and obtain the soft decisions (line 1), i.e., for ever i ∈ [n] we get the 4 numbers (p i,i , p i,x , p i,y , p i,z ), where p i,E can be interpreted as the probability P(e i = E) of the error E ∈ {i, x, y, z} in the i-th qubit. Next, we find the hard decisions for all qubits by the formula (line 2): If the BP decoding is successful (i.e., we have the correct syndrome vector), then there is no need in the postprocessing, and the BP-OSD decoder returns the result of the BP decoder (line 4).
After the first stage, the probabilities p i,E are also used in the BP-OSD algorithm to sort the qubits in the order of increasing reliability ρ(p i,i , p i,x , p i,y , p i,z ), where we can define the reliability of the i-th qubit as However, in all our simulations we used the formula which gives similar error-correcting performance, but it is a little bit easier to calculate. After we sorted the qubit positions (line 7), the qOSD algorithm from the previous subsection is used to recover the errors in the least reliable qubits from the hard decisions for the remaining qubits (line 8).

end
In Algorithm 3, we used the qOSD algorithm, which is a good choice for the standard depolarizing channel. If we have a channel with independent X and Z errors, then for a CSS code defined by a pair of parity-check matrices H X , H Z we can also use the standard OSD (Algorithms 1 and 2) for the X and Z components of the error vector e separately. In this case, lines 6-8 in Algorithm 3 can be replaced by the following steps: 1. Calculate the X and Z error probabilities: 2. Sort the qubits in the decreasing order of their X error and Z error probabilities and let σ X , σ Z ∈ S n be the corresponding permutations such that: Run the OSD decoder separately for the X and Z components ofê:  For all our simulations we used the normalized minsum (NMS) decoder with the normalization factor 0.625, which approximates the non-binary BP decoder from [21], [11,Algorithm 1] in log-domain and is more numerically stable in some cases. The maximal number of iterations was set to 32. We used the layered scheduling in order to increase the convergence speed of the decoder by approximately two times 11 . For a good review of practical aspects of the BP decoder implementation, see [35,Chapter 4]. The error-correcting performance in our simulations is measured either in terms of the logical error rate or the word error rate 12 (WER). We should also stress that in many cases we use only the OSD-0 algorithm (Algorithm 2), which has complexity O(n 3 ), while the complexity in the general case is O(n 3 + n2 w ).
In Fig. 1 we show the effect of the OSD-0 postprocessing after the BP decoder for different QLDPC codes, including some known and new codes described in the next sections. As you can see, the gain of the BP-OSD over the BP decoder in terms of WER for some codes is up to 5 orders of magnitude (code B1). However, for code A1 the difference between BP-OSD and BP is quite small. We think that this is because the column weight of the matrices H X and H Z is quite high (it is equal to 5), and hence the performance of the BP decoder is very good even without the post-processing.
Remark. It is interesting to note that 8-limited Haah's [[1024,30]] code, which has local stabilizers in 3D, also performs very well under the BP-OSD decoder. To the best of our knowledge, this is the first such demonstration on the depolarizing channel. From our point of view, this observation implicitly suggests that some of Haah's cubic codes may have very good minimum distances. In fact, even after very long runs of the BP-OSD decoder on this [[1024,30]] code, we have not found any non-degenerate codewords of weight less than 32.
In Fig. 2 we show the impact of the OSD order w on the WER performance of the 8-limited [[882,48,16]] code (Appendix B, code B2) under the BP-OSD decoder. We see that the WER performance in this case can be further improved by increasing the OSD order w. At the same time, from our experiments with the OSD postprocessing on different QLDPC codes, we observed that in many cases the impact of the OSD order on the WER performance is quite small.

Different post-processing algorithms
In this section, the OSD post-processing algorithm is compared against the other known post-processing methods that also significantly improve the performance of the BP decoder. We assume that the prob-ability distributions in the BP decoder are represented in terms of the log-likelihood ratios (LLRs). Below you can find a brief description of these methods, where each method can be applied repeatedly after the BP decoder until all the parity-checks are satisfied or the maximal number of attempts n a is reached.
• Random perturbation [21]. If the syndrome is non-zero after the BP decoding, then we randomly choose an unsatisfied c-node and randomly change the initial LLRs for all the v-nodes adjacent to this c-node. The main parameter of this method is the variance of the perturbation magnitude. After the changes are done, the BP decoder runs with the perturbed input LLRs.
• Enhanced feedback [22]. It is similar to the previous approach but the perturbations are not random and calculated using the previous BP output. If after the BP decoding the parity-checks are not satisfied, we randomly choose an unsatisfied c-node. Then for all the v-nodes adjacent to this cnode, we set the initial LLRs using the BP decoder output for these nodes. After this is done the BP decoder runs with the perturbed input LLRs.
• Matrix augmentation [23]. In this method instead of modification of the input LLRs, the paritycheck matrix itself is modified by random duplication of some rows. The fraction δ of the duplicated rows is called the augmentation density. Then a new decoding attempt with the augmented matrix is performed.
To compare the performance of all the described methods with the BP-OSD decoder we use the 6-limited [[1270, 28]] QLDPC code (B3 in Appendix B). This code belongs to the new class of codes described in Section 5. For all these methods we used n a = 100. We can see in Fig. 3 that the OSD post-processing outperforms all the above-mentioned post-processing methods and also outperforms the 4-limited [[1201, 1,25]] surface code on the MPS decoder from [5], which is almost optimal for this code. Let us note that all the other post-processing methods also have the WER gain that is more than 10 3 over the BP decoder for code B3. In fact, we observed a similar WER gain for many other codes from the class of (3, 6)-regular CSS QLDPC codes with a sufficiently large minimum distance. We think that this is mainly because the classical (3, 6)-regular LDPC codes defined by the matrices H X and H Z of such CSS codes themselves have many harmful trapping sets. Thus these (3, 6)-regular QLDPC codes additionally have a lot of degenerate codewords of very low weight starting from 6.

New generalized bicycle codes 4.1 Ansatz with commuting matrices
The commutativity conditions such as (4) and (5) are a serious obstacle to designing good QLDPC codes using random-like constructions similar to the constructions used for classical LDPC codes. Thus it makes sense to consider large families of matrices of some particular form called ansatz, where the commutativity conditions are always satisfied. One such quite general ansatz for CSS codes was proposed in [10] as a generalization of the bicycle QLDPC codes [37]. Let us briefly remind this ansatz. Consider two commuting binary n × n matrices A and B, i.e., AB = BA. Let Then we see that H X H T Z = AB + BA = 0, and the commutativity condition (5) is always satisfied. It was proposed in [10] to use binary circulant matrices A and B since they always commute. The corresponding class of codes is called the generalized bicycle (GB) codes, where the bicycle codes [37] are obtained as a special case when B = A T . where a 0 , . . . , a −1 ∈ F q . It is readily seen that the matrix A can be represented in the form

Ring of circulants
where I is the × identity matrix and is the × permutation matrix representing the right cyclic shift by one position. Since P = I, we see that the ring of all × circulant matrices over F q is isomorphic to the ring F q = F q [x]/(x − 1) of polynomials over F q modulo the polynomial x − 1.
Hence the circulant matrix A can be uniquely represented by the polynomial a(x) = a 0 + a 1 x + · · · + a −1 x −1 and the product C = AB of two circulant matrices represented by polynomials a(x), b(x) ∈ R corresponds to the polynomial

Dimension of generalized bicycle codes
As we saw before, to define two binary circulant × matrices A and B we need to provide two binary poly-  (x), b(x) ∈ F 2 is given by the formula: where g(x) = gcd(a(x), b(x), x − 1). 13 Let us point out that this dimension formula was given in the paper [10, Theorem 2] in a slightly more complex form. A similar formula was proved only in the special case of single generator codes.
To prove this formula we show that rk H X = rk H Z = n − deg g(x) and use CSS dimension formula (6).
Let us first find the rank of the matrix H X = [A, B], which is equal to the dimension of its column space. It is easy to see that the column space of H X (called its syndrome space) is equal to the following set: Using the described above polynomial representation of column vectors and circulant matrices we can consider this set as the following set of polynomials from F 2 : It is easy to verify that this set is the principal 14 ideal of the ring F 2 generated by g(x). Hence it is the cyclic code C g with generator polynomial g(x), and we proved that rk H X = dim C g = n − deg g(x). We call this cyclic code C g the syndrome code of H X since its codewords are precisely the polynomial representations of the syndrome space of H X .
To complete the proof we also need to show that rk H Z = n − deg g(x). Using similar arguments as above we can consider the syndrome code for H Z and show that it is generated by the "transposed" polynomial g * (x) = g(x −1 ). Though the codes generated by the polynomials g(x) and g * (x) may differ, they always have the same dimension since the corresponding circulant matrices G and G T have the same rank. Hence we also proved that rk H Z = n − deg g(x), and the proof of formula (14) is complete.

Construction methods
The proof of Proposition 1 gives us also some valuable information on how to find generalized bicycle codes of high dimension. If we fix the circulant size then all possible dimensions k of the generalized bicycle codes with this circulant size are characterized by the degrees of all possible factors of the polynomial x − 1. Indeed, for each factor g(x) of the polynomial x − 1 we can always choose polynomials a(x), b(x) ∈ F 2 such that: since these polynomials are just two codewords from the cyclic code C g generated by g(x), which we called the syndrome code of H X . To produce a w-limited QLDPC we just need to find low weight polynomials a(x) and b(x) that are the codewords of C g . In practice, this can be accomplished by several methods. If the circulant size is relatively small, we can find a(x), b(x) by an exhaustive search over all polynomials of the given weight.
Another alternative is to generate random polynomials of a given weight from F 2 until we find a pair of polynomials that satisfies condition (15). Since the probability that a random polynomial of a given weight belongs to C g is equal approximately to 2 − deg g (x) ; then if we test more then 2 deg g(x) random polynomials we will find a polynomial from C g with high probability.
When we find a pair of polynomials a(x), b(x) we also need to check that the corresponding code has good error correcting performance. This can be done by a simulation of the corresponding generalized bicycle code.
All the generalized bicycle codes from Appendix B were found by the described above methods.

GB codes with syndrome protection
Another important observation, made in the proof of Proposition 1, is that the syndrome codes of the paritycheck matrices H X and H Z are the cyclic codes with the generator polynomials g(x) and g * (x), respectively. Let us mention that the syndrome code of a paritycheck matrix is precisely the set of all possible syndromes for it. Hence if we use generator polynomials g(x) and g * (x) that define cyclic codes with minimum distance d, we see that the syndromes for matrices H X and H Z are protected by these cyclic codes. Since the syndrome measurements for quantum codes are performed by faulty hardware, some additional protection of the syndromes may be used to improve the reliability of the syndrome measurements [38,39]. Let us also mention that the polynomials g(x) and g * (x) always produce cyclic codes C g and C g * with the same minimum distance since the "transpose" map is an automorphism of the ring F 2 that respects the weight of the polynomials, and therefore we have that Hence we can use any cyclic code with generator polynomial g(x) and minimum distance d to construct a CSS code, where the syndromes for H X and H Z are protected by cyclic codes of minimum distance d. It is important to note that since a(x), b(x) ∈ C g , the weight of the polynomials a(x), b(x) can not be smaller than this minimum distance d.

Comparison with other codes
In this subsection, we consider several new examples of generalized bicycle codes and compare their perfor-  1,9]] homological product (HMP) code from [8] under the BP decoder with the OSD-10 post-processing.
mance under the BP-OSD decoder against some other already known QLDPC codes.
Example 1. Let us consider the primitive narrow-sense BCH [127, 14,5] code C g with the generator polynomial: If we set a(x) = 1 + x 15 + x 20 + x 28 + x 66 and b(x) = 1 + x 58 +x 59 +x 100 +x 121 , then we obtain the 10-limited generalized bicycle [[254, 28]] code. Its minimum distance is not available, but the performance of this code (see Fig. 1 The definition of these codes is given in Section 5. If we set a(x) = 1 + x + x 14 + x 16 + x 22 , b(x) = 1 + x 3 + x 13 + x 20 + x 42 we obtain the 10-limited generalized bicycle [[126, 28]] code. Its performance with the standard BP decoder (binary and non-binary) is shown in Fig. 5, code A2. We compared its performance with the performance of the neural BP decoder for the bicycle [[256, 32]] code constructed in [41]. Such a big difference in the performance is mostly because the QLDPC [[256, 32]] code used in [41] has a small minimum distance compared to the [[126, 28]] code, which minimum distance is 8. This minimum distance was found by an exhaustive search (see Appendix B). Another possible reason is that the neural BP decoder proposed in [41] was based on the binary BP, which usually has worse performance than its non-binary version.  code A3 and the [[46, 2, 9]] code A4, see Appendix B). We compared their performance (see Fig. 4) with the performance of an 8-limited [[49, 1, 9]] homological product (HMP) code from [8] using the BP with OSD-like post-processing. As we can see the performance of the newly constructed codes is similar to the [[49, 1,9]] code. At the same time, their rates are higher.
Example 4. In Fig. 6 14]] code. At the same time, it has the same weight of stabilizers, and its code length is 5 times smaller.

Hypergraph product (HP) codes
In this section, we propose a new generalization of hypergraph product codes [7] in the case when one of the parity-check matrices in the product is square. Let us first remind the definition of these codes in a matrix form [9]. Suppose we have an [n a , k a , d a ] linear code C a and an [n b , k b , d b ] linear code C b defined by parity-check matrices 17 a ∈ M ma×na (F 2 ) and b ∈ M m b ×n b (F 2 ) respectively. Then the hypergraph product code is the −m a ). As it was shown in [7], the minimum distance d of the hypergraph product code C satisfies the following lower bound: are the minimal distances of the "transposed" codes C T a and C T b defined by the paritycheck matrices a T and b T respectively. It is important to note that if the matrices a and b are w-limited, then the corresponding CSS code C is 2w-limited. Hence, using known asymptotically good families of classical LDPC codes with (w c , w r )-limited parity checkmatrices, it is possible [7] to construct w-limited CSS codes with asymptotically non-zero rate and d = Θ( √ n) as n → ∞. In [9] the hypergraph product construction was further improved, and it was shown that one can construct good hypergraph product codes using square parity-check matrices a and b. In fact, many of the best-known small-length hypergraph product codes are constructed using square parity check matrices of cyclic codes (see [9,10]). In [10] hyperbicycle CSS codes, which generalize both generalized bicycle and hypergraph product codes, were proposed. Here we consider another generalization of hypergraph product codes where the matrix b is square.

Generalized hypergraph product codes
In what follows by a ring we always mean a ring with identity. Let R be a ring. We denote the ring of all m × n matrices over R by M m×n (R) or by M n (R) when m = n. If R is the ring of × matrices over some field F we identify the elements of M m×n (R) with the corresponding block matrices from M m ×n (F).
Consider a binary matrix b ∈ M (F 2 ). We say that the matrix b and a ring R ⊆ M (F 2 ) commute if all matrices from R commute with b. Example 5. Consider b ∈ M (F 2 ) and R = {0, I}, where 0 and I are the zero and the identity matrices from M (F 2 ) respectively; then b and R always commute. Example 6. Let b be a binary circulant matrix and R be the ring of all binary circulant matrices of the same size; then b and R always commute.
Suppose that a matrix b ∈ M (F 2 ) and a ring R ⊆ M (F 2 ) commute. Consider a matrix A = (a ij ) m×n ∈ M m×n (R). We denote by C (A, b) the CSS code called a generalized hypergraph product (GHP ) code with the following parity-check matrices 18 : 18 Let us warn the reader that we understand H X and H Z as where A * = (a T ji ) n×m , and I m , I n are the identity matrices over R of size m and n respectively. The correctness of this definition follows from the following: The code length N of the CSS code C (A, b) is given by N = (m + n) . We will show later how to find the dimension K of the code C (A, b) in a special case.
One can easily verify that if we take a matrix b ∈ M (F 2 ) and the ring R = {0, I} (as in Example 5) then the CSS code C (A, b) is the hypergraph product code defined by the binary m × n matrixã = (ã ij ) m×n and the binary × matrix b, where for all i ∈ [m], j ∈ [n] we have:ã We can also see that the ansatz with two commuting matrices A and B given in (12) is also a special case of the new ansatz described in (16), where the matrix A is a 1 × 1 block matrix.

Quasi-cyclic generalized hypergraph product codes
In this subsection, we describe quasi-cyclic (QC) GHP codes. In fact, it can be shown that this particular subclass of GHP codes is equivalent (up to some permutation of qubits) to a special case of hyperbicycle codes from [10,Eq. (19), χ = r 2 = n 2 = 1]. However, we believe that the concise language of polynomial matrices adopted in the current work is much more suitable for our examples, which are essentially sparse random QC GHP codes subject to some additional constrains. Now, let us take b and R as in Example 6. Hence b is a binary circulant matrix, and R is the ring of all binary circulant matrices of the same size. In this case, the matrices H X and H Z defined by (16) are block matrices, where each block is a binary circulant matrix of size . Such matrices are called quasi-cyclic. Let us note that quasi-cyclic (QC) matrices are well known in classical coding theory. In fact, most of the best-known practical classical LDPC codes have QC parity-check matrices. We will show in this section how to find the dimension of generalized hypergraph product codes defined by (16). For simplicity, we consider here only the case when the circulant size is odd. Before we can provide the formula for the dimension we need some supplementary definitions from algebra. the corresponding binary block matrices (not as matrices over M (F 2 )).
Here we adopt the polynomial representation of QC matrices used in [42,43]. For any polynomial p(x) ∈ F q [x] of degree d we consider the ring with addition and multiplication modulo p(x). By F d q we denote the d-dimensional space of the d × 1 column vectors over F q . We also identify an element f (x) ∈ F q [x]/(p(x)) with the corresponding column vector f = (f 0 , . . . , f d−1 ) ∈ F d q . Let us recall that by F 2 we denote the ring of circulants F 2 [x]/(x − 1). We use the standard identification of the circulant × matrices over F 2 with the elements of the ring F 2 (see Subsection 4.2), where a column vector a ∈ F 2 corresponds to the circulant matrix with the first column equal to a. Using this identification we can consider an m × n QC matrix over F 2 of circulant size as an m × n matrix over the ring F 2 . We also consider n × 1 column vectors over F 2 as n × 1 column vectors over F 2 . Given the above identification we consider multiplication of an m × n QC matrix by an n × 1 column vector over F 2 as multiplication of an m × n matrix by an n × 1 column vector over F 2 .
The algebraic structure of the ring F 2 is very well studied in the literature (see Appendix A for further details). Since we consider the case when is odd, the polynomial x − 1 factors into a product of irreducible polynomials over F 2 : Hence the ring F 2 is isomorphic to the direct product of finite fields: where the field Let us consider the maps ϕ i : F 2 → F i given by the formula: We also naturally extend this map to any vectors and matrices over F 2 . The following lemma is the key to the dimension formula for the QC generalized hypergraph product codes.

Lemma 1.
Let A ∈ M m×n (F 2 ) be a binary QC matrix. Then its rank (over F 2 ) is given by: Proof. The lemma easily follows from the isomorphism shown in (18) between F 2 and the direct product of where u 1 , . . . , u n ∈ F 2 . Therefore we see that the cardinality of the column space of the non-binary matrix ϕ i (A) over the field F i is equal to (2 di ) ri = 2 diri . Hence using isomorphism (18) we conclude that the number of different vectors in the column space of the matrix A is equal to and we proved that The following proposition provides the formula for the dimension of the QC CSS code C(A, b) if is odd. since for all i ∈ [s] \ S we have that ϕ i (b) = 0 and therefore the non-binary matrix ϕ i (H X ) is full rank. Hence by applying Lemma 1 to the matrix H X we have: To complete the proof we need also to find the rank of the matrix H Z . It is easier to find the rank of the matrix H Z obtained from H Z by the application of the "transpose" map u → u * to each element. Since the transpose map is an automorphism on the ring F 2 we see that the number of vectors in the row space of H Z and the row space of H Z is the same. Hence rk H Z = rk H Z . The rank of the matrix H Z = [bI n , A T ], can be found in the same way as for the matrix H X : Now we apply formula (6) for the CSS dimension and obtain: This concludes the proof.
Let us note that if the polynomial b(x) is such that g(x) = gcd(b(x), x −1) is an irreducible factor of x −1, then we can give a more elegant formula for the dimension of the code C(A, b) in some special cases. Indeed, in this case we have: where F is the finite field F 2 [x]/(g(x)) and ϕ(A) is the F -image of A under the action of the map ϕ : u(x) → u(x) mod g(x). Since ϕ(A) is a matrix over F , it defines, as a parity-check matrix, a non-binary linear code over F of dimension At the same time, b(x) defines the cyclic code as a check polynomial. Thus formula (19) gives us the dimension K of C (A, b) in terms of the dimensions k A and k b for the two special cases shown below: 1. if A is a square matrix, i.e. m = n, we have In the current work, we construct only codes that correspond to case 1 above. Such codes are defined by an irreducible factor b(x) of x − 1 and a square matrix A over F 2 . To construct all the examples of QC GHP codes given in Appendix B (codes B1, B2, and B3), we used the following procedure. First we fix a low weight irreducible polynomial b(x) ∈ F 2 [x] such that b(x) | x − 1, and then randomly choose a polynomial matrix A subject to some constraints. Since we want to obtain sparse parity-check matrices H X and H Z for the code C (A, b), we restrict each entry of A to be either 0 or a monomial x i , 0 i < . When the number of non-zero entries in each row and each column of A is bounded above by w, this restriction guaranties that both H X and H Z are (w + deg b(x))-limited matrices. Another restriction is related to the girth of the Tanner graphs T X , T Z for H X , H Z . As you can see from Tab. 1 in Appendix B, the girth of the Tanner graphs T X , T Z for all the constructed QC GHP codes is equal to 6, which means that T X , T Z do not contain 4-cycles. This restriction helps to improve the error-correcting performance under the BP decoder.
Example 7. In Fig. 7 you can see the WER performance (under the BP-OSD decoder) of the two codes: the 6-limited [[1922,50,16]] hypergraph product (HP) code (C2 in Appendix B) and the 6-limited [[882,24]] generalized hypergraph product (GHP) code (B1 in Appendix B). You can see from Fig. 7 that the HP [[1922,50,16]] code has some error floor even under the BP-OSD decoder. We believe that this is due to a large amount of low-weight non-degenerate codewords in this code.

Conclusion
We proposed new OSD-like post-processing for the BP decoder that shows on some codes much better performance than all the modifications known to the authors. We also constructed several new generalized bicycle codes that show very good performance compared to the other known codes with similar parameters. We proposed a new ansatz for quantum CSS codes and showed how to estimate the dimension of such codes in some special cases. Unfortunately, we have not found any nontrivial general lower bound on the minimum distance of such codes. We think that to find such a bound is an interesting open problem 19 since this class contains 19 Since the first version of this work appeared in arXiv, some lower bounds on the minimal distance have been obtained in [16,17,18] for special cases of QC GHP codes C(A, 1 + x). some of the best known QLDPC codes, and their practical performance under the BP-OSD decoder is also quite good. Finally, we compared the performance of one of our codes from the new family and showed that it has better performance than a relatively large surface code of similar code length even if this code is decoded by a near-optimal decoder.

B Matrices used for simulations
All matrices for generalized hypergraph product and generalized bicycle codes that we used for simulations have the form H X = [A, B], H Z = [B T , A T ] where A and B are quasi-cyclic matrices. Thus, to define the code, we specify matrices A and B in the polynomial form as matrices over F 2 .
For all codes presented here, we also provide lower and upper bounds on the minimum distance obtained either by methods similar to the ones from [45] or by a straightforward reduction of the minimum distance problem to a mixed integer linear program and using the GNU Linear Programming Kit, Version 4.63, http://www.gnu.org/software/glpk/glpk.html.  C. Hypergraph product (HP) codes. Each hypergraph product code in our simulations is constructed from a single cyclic code defined by its parity polynomial h(x) and the length .

C Additional Simulations
In this appendix, we show some additional simulation results of the BP-OSD decoder on the depolarizing channel. In Fig. 9, we demonstrate the WER performance of several 6-limited generalized bicycle (GB) codes with the parameters [[2 s+1 − 2, 2s]], s ∈ N, under the BP-OSD-0 decoder. As one can see from the curves, these codes have the threshold p GB ≈ 15%, which is quite close to the corresponding threshold p S ≈ 18% of the surface codes from [5]. The codes from this family were constructed in the same way as the GB codes from Subsection 4.6, where for each s ∈ N the corresponding syndrome code C g is the cyclic Hamming [2 s − 1, s, 3] code defined by some irreducible polynomial g(x) ∈ F 2 [x] of degree s.
In all our previous examples, we considered only CSS codes. In fact, as it is shown in Subsection 3.3, the BP-OSD decoder can be applied to any stabilizer code. In Fig. 10, we demonstrate the WER performance of the BP-OSD-0 decoder on the 5-limited cyclic non-CSS [[126,2,12]] code. This code is defined by the paritycheck matrix H = [H X | H Z ], where H X and H Z are the 126×126 circulant matrices represented respectively by the polynomials: 1 + x 71 + x 55 and 1 + x 40 + x 86 . As we can see, the WER gain of OSD-0 post-processing in this case is at least three orders of magnitude.  Figure 9: The WER performance on the depolarizing channel for five 6-limited GB codes under the BP-OSD-0 decoder. The threshold is p GB ≈ 15%. The corresponding threshold of the 4-limited surface codes [5] is p S ≈ 18%.