# Introduction Low density parity check (LDPC) codes are among the best performing codes in practice, being capable of correcting errors close to the Shannon limit. Magma provides facilities for the construction, decoding, simulation and analysis of LDPC codes. ## Constructing LDPC Codes LDPC codes come in two main varieties, *regular* and *irregular*, defined by the row and column weights of the sparse parity check matrix. If all columns in the parity check matrix have some constant weight $a$, and all rows have some constant weight $b$, then the LDPC code is said to be $(a,b)$-regular. When either the columns or the rows have a distribution of weights, the LDPC code is said to be irregular. Currently, there do not exist many techniques for the explicit construction of LDPC codes. More commonly, these codes are selected at random from an ensemble, and their properties determined through simulation. ### `LDPCCode(H): MtrxSprs -> Code` Given a sparse binary matrix $H$, return the LDPC code which has $H$ as its parity check matrix. ### `GallagerCode(n, a, b): RngIntElt, RngIntElt, RngIntElt -> Code` Return a random $(a,b)$-regular LDPC code of length $n$, using Gallager’s original method of construction. The row weight $a$ must divide the length $n$. ### `RegularLDPCEnsemble(n, a, b): RngIntElt, RngIntElt, RngIntElt -> Code` Return a random code from the ensemble of $(a,b)$-regular binary LDPC codes. ### `IrregularLDPCEnsemble(n, Sc, Sv): RngIntElt, SeqEnum, SeqEnum -> Code` Given (unnormalized) distributions for the variable and check weights, return length $n$ irregular LDPC codes whose degree distributions match the given distribution. The arguments $Sv$ and $Sc$ are sequences of real numbers, where the $i$-th entry indicates what percentage of the variable (resp. check) nodes should have weight $i$. Note that the distributions will not be matched perfectly unless everything is in complete balance. ### `MargulisCode(p): RngIntElt -> Code` Return the (3,6)-regular binary LDPC code of length $2(p^3-p)$ using the group-based construction of Margulis. ### `Example: Is LDPC (ex-f1c0a4)` Most LDPC constructions are generated pseudo-randomly from an ensemble, so the same function will return a different code each time. To be able to re-use an LDPC code, the sparse parity check matrix must be saved. ```magma > C1 := RegularLDPCEnsemble(10, 2, 4); > C2 := RegularLDPCEnsemble(10, 2, 4); > C1 eq C2; false > LDPCMatrix(C1):Magma; SparseMatrix(GF(2), 5, 10, [ 4, 2,1, 3,1, 4,1, 6,1, 4, 1,1, 7,1, 9,1, 10,1, 4, 1,1, 2,1, 3,1, 7,1, 4, 5,1, 6,1, 8,1, 9,1, 4, 4,1, 5,1, 8,1, 10,1 ]) > H := LDPCMatrix(C1); > C3 := LDPCCode(H); > C3 eq C1; true ``` ## Access Functions Since a code can have many different parity check matrices, the matrix which defines a code as being LDPC must be assigned specifically. Any parity check matrix can be assigned for this purpose, and once a code is assigned an LDPC matrix it is considered by Magma to be an LDPC code (regardless of the density or other properties of the matrix). The matrix must be of sparse type (`MtrxSprs`). ### `IsLDPC(C): Code -> BoolElt` Return true if $C$ is an LDPC code (which is true if it has been assigned an LDPC matrix). ### `AssignLDPCMatrix(~C, H): Code, MtrxSprs` Given a sparse matrix $H$ which is a parity check matrix of the code $C$, assign $H$ as the LDPC matrix of $C$. ### `LDPCMatrix(C): Code -> MtrxSprs` Given an LDPC code $C$, return the sparse matrix which has been assigned as its low density parity check matrix. ### `LDPCDensity(C): Code -> FldReElt` Given an LDPC code $C$, return the density of the sparse matrix which has been assigned as its low density parity check matrix. ### `IsRegularLDPC(C): Code -> BoolElt` Returns `true` if $C$ is an LDPC code and has regular column and row weights. If true, the row and column weights are also returned. ### `TannerGraph(C): Code -> Grph` For an LDPC code $C$, return its Tanner graph. If there are $n$ variables and $m$ checks, then the graph has $n+m$ nodes, the first $n$ of which are the variable nodes. ### `LDPCGirth(C): Code -> RngIntElt` For an LDPC code $C$, return the girth of its Tanner graph. ### `LDPCEnsembleRate(v, c): RngIntElt, RngIntElt -> FldReElt` ### `LDPCEnsembleRate(Sv, Sc): [FldReElt], [FldReElt] -> FldReElt` Return the theoretical rate of LDPC codes from the ensemble described by the given inputs. ### `Example: Is LDPC (ex-69eb66)` In Magma, whether or not a code is considered LDPC is based solely on whether or not an LDPC matrix has been assigned. This example shows that any code can be made to be considered LDPC, although a random parity check matrix without low density will perform very badly using LDPC decoding. ```magma > C := RandomLinearCode(GF(2),100,50); > IsLDPC(C); false > H := SparseMatrix(ParityCheckMatrix(C)); > H; Sparse matrix with 50 rows and 100 columns over GF(2) > AssignLDPCMatrix(~C, H); > IsLDPC(C); true > LDPCDensity(C); 0.253400000000000014122036873232 ``` The density of the parity check matrices of LDPC codes is much lower than that of randomly generated codes. ```magma > C1 := RegularLDPCEnsemble(100,3,6); > C1:Minimal; [100, 50] Linear Code over GF(2) > LDPCDensity(C1); 0.0599999999999999977795539507497 ``` ## LDPC Decoding and Simulation The impressive performance of LDPC codes lies in their iterative decoding algorithm. Magma provides facilities to decode using LDPC codes, as well as simulating transmission over a binary symmetric or white Gaussian noise channels. The binary symmetric channel transmits binary values and is defined by $p < 0.5$. Each individual bit independently sustains a “bit-flip” error with probability $p$. The Gaussian channel is analog, transmitting real-values, and is defined by a standard deviation $\sigma$. Binary values are mapped to $-1$ and $1$ before being transmitted. Each value independently sustains an errors which are normally distributed about $0$ with standard deviation $\sigma$. ### `LDPCDecode(C, v): Code, ModTupRngElt -> ModTupRngElt` ```magma Channel : MonStgElt Default: "BinarySymmetric" p : RngReSubElt Default: 0.1 StdDev : RngReSubElt Default: 0.25 Iterations: RngIntElt Default: Dimension(C) ``` For an LDPC code $C$ and a received vector $v$, decode $v$ to a codeword of $C$ using the LDPC iterative decoding algorithm. The nature of the channel from which $v$ is received is described by the variable argument `Channel`, which can either be the `BinarySymmetric` channel or the `Gaussian` channel. Errors on the binary symmetric channel is described by the argument `p`, while on the Gaussian channel they are described by `StdDev`. The vector $v$ must be over a ring corresponding to the channel which is selected. For the binary symmetric channel $v$ must be a binary vector over $F_2$, while for the Gaussian channel it must be real-valued. Since the decoding algorithm is iterative and does not necessarily terminate on its own, a maximum number of iterations needs to be specified using the argument `Iterations`. The default value is much larger than would normally be used in practice, giving maximum error-correcting performance (at possibly some cost to efficiency). ### `Example: DecodeLDPC BSC (ex-ece197)` Errors in the binary symmetric channel are just bit flips. ```magma > n := 500; > C := RegularLDPCEnsemble(n, 4, 8); > e := 5; > Errs := {}; > repeat Include(~Errs, Random(1,n)); until #Errs eq e; > v := Random(C); > ev := AmbientSpace(C)![ (i in Errs) select 1 else 0 : i in [1..n]]; > rec_vec := v + ev; > time res := LDPCDecode(C, rec_vec : Channel:="BinarySymmetric", p:=0.2); Time: 0.000 > res eq v; true ``` ### `Example: DecodeLDPC BSC (ex-557473)` For the Gaussian channel binary vectors are considered to be transmitted as sequences of the values $1$ and $-1$. Errors are normally distributed with a standard deviation defined by the channel. To simulate a Gaussian channel requires obtaining normally distributed errors. This can be done (discretely) by generating a multiset of possible errors. ```magma > sigma := 0.5; > MaxE := 3.0; > N := 100; > V := [ MaxE*(i/N) : i in [-N div 2..N div 2]]; > E := [ 0.5*(1+Erf(x/(sigma*Sqrt(2)))) : x in V ]; > Dist := {* V[i]^^Round(1000*(E[i]-E[i-1])) : i in [2..#V]*}; ``` A codeword of an LDPC code needs to be mapped into the real domain. ```magma > n := 500; > C := RegularLDPCEnsemble(n, 4, 8); > v := Random(C); > R := RealField(); > RS := RSpace(R, n); > vR := RS ! [ IsOne(v[i]) select 1 else -1 : i in [1..n]]; ``` Normally distributed errors are then introduced, and the received vector decoded. ```magma > for i in [1..n] do > vR[i] +:= Random(Dist); > end for; > time res := LDPCDecode(C, vR : Channel:="Gaussian", StdDev:=sigma); Time: 0.000 > res eq v; true ``` ### `LDPCSimulate(C, N): Code, RngIntElt -> FldReElt, FldReElt` ```magma Channel : MonStgElt Default: "BinarySymmetric" p : RngReSubElt Default: 0.1 StdDev : RngReSubElt Default: 0.25 Iterations: RngIntElt Default: Dimension(C) ``` For an LDPC code $C$, simulate $N$ transmissions across the given channel and return the accumulated bit error rate and word error rate. The variable arguments are as described for the function `LDPCDecode`. The channel which is described controls not only the way the decoding algorithm functions, but also the nature of the errors introduced during the simulation. ### `Example: Decode Ensemble (ex-103564)` The more noise that is introduced into the channel the error rate increases. Note that the bit error rate (the first return value) is always lower than the word error rate (the second return value). ```magma > C := RegularLDPCEnsemble(200, 3, 6); > LDPCSimulate(C, 10000 : Channel := "Gaussian", StdDev := 0.7); 0.00118249999999999995739519143001 0.00619999999999999978211873141731 > LDPCSimulate(C, 10000 : Channel := "Gaussian", StdDev := 0.75); 0.00749499999999999992617016886243 0.0410000000000000017208456881690 > LDPCSimulate(C, 10000 : Channel := "Gaussian", StdDev := 0.8); 0.0337220000000000019735324485737 0.159700000000000008615330671091 > LDPCSimulate(C, 10000 : Channel := "Gaussian", StdDev := 0.85); 0.0856174999999999991606713933834 0.370800000000000018474111129763 > LDPCSimulate(C, 10000 : Channel := "Gaussian", StdDev := 0.9); 0.162790999999999991265653420669 0.640499999999999958255614274094 > LDPCSimulate(C, 10000 : Channel := "Gaussian", StdDev := 0.95); 0.237657499999999993756105709508 0.840099999999999957900342906214 > LDPCSimulate(C, 10000 : Channel := "Gaussian", StdDev := 1.0); 0.296526500000000026169288958044 0.944799999999999973177011725056 ``` ## Density Evolution The asymptotic performance of ensembles of LDPC codes can be determined using *density evolution*. An ensemble of LDPC codes (either regular or irregular) is defined by a pair of degree distributions, corresponding to the degrees at the variable and check nodes of the Tanner graph. Over a specific channel, the critical parameter which defines the asymptotic performance of a given ensemble is its *threshold*, which is a value of the channel parameter (i.e., the probability of error $p$ for the binary symmetric channel, and the standard deviation $\sigma$ for the Gaussian channel). When the channel parameter is less than the threshold, asymptotically a code from the ensemble will decode with an error probability tending to zero. However, at any channel parameter above the threshold there will be a non-vanishing finite error probability. Determining the threshold of an ensemble over the binary symmetric channel is relatively trivial, however over the real-valued Gaussian channel it can involve extensive computations. The speed depends heavily on the granularity of the discretization which is used, though this also affects the accuracy of the result. The default settings of the Magma implementation use a reasonably coarse discretization, emphasizing speed over accuracy. These (still quite accurate) approximate results can then be used to help reduce the workload of calculations over finer discretizations if more accuracy is required. ### `LDPCBinarySymmetricThreshold(v, c): RngIntElt, RngIntElt -> FldReElt` ### `LDPCBinarySymmetricThreshold(Sv, Sc): [FldReElt], [FldReElt] -> FldReElt` ```magma Precision: RngReSubElt Default: 0.00005 ``` Determines the threshold of the described ensemble of LDPC codes over the binary symmetric channel, which is the critical value of the channel parameter above which there is a non-vanishing error probability (asymptotically). The ensemble can either be defined by two integers for $(v,c)$-regular LDPC codes, or by two density distributions $Sv$ and $Sc$, which are sequences of non-negative real numbers. The density distributions do not need to be normalized, though the first entry (corresponding to weight $1$ nodes in the Tanner graph) should always be zero. The computation proceeds by establishing lower and upper bounds on the threshold, then narrowing this range by repeatedly performing density evolution on the midpoint. The argument `Precision` controls the precision to which the threshold is desired. ### `DensityEvolutionBinarySymmetric(v, c, p): RngIntElt, RngIntElt, FldReElt -> BoolElt` ### `DensityEvolutionBinarySymmetric(Sv, Sc, p): [FldReElt], [FldReElt], FldReElt -> BoolElt` Perform density evolution on the binary symmetric channel using channel parameter $p$ and determine the asymptotic behaviour for the given LDPC ensemble. The return value is boolean, where `true` indicates that $p$ is below the threshold and the ensemble has error probability asymptotically tending to zero. ### `Example: DE BSC (ex-b0bbfa)` Density evolution on the binary symmetric channel is not computationally intensive. ```magma > time LDPCBinarySymmetricThreshold(3, 6); 0.0394140625000000000000000000000 Time: 0.010 > time LDPCBinarySymmetricThreshold(4, 8); 0.0473925781250000000000000000001 Time: 0.110 > time LDPCBinarySymmetricThreshold(4, 10); 0.0368359375000000000000000000000 Time: 0.090 ``` ### `LDPCGaussianThreshold(v, c): RngIntElt, RngIntElt -> FldReElt` ### `LDPCGaussianThreshold(Sv, Sc): [FldReElt], [FldReElt] -> FldReElt` ```magma Lower : RngReSubElt Default: 0 Upper : RngReSubElt Default: Infinity() Points : RngIntElt Default: 500 MaxLLR : RngReSubElt Default: 25 MaxIterations: RngIntElt Default: Infinity() QuickCheck : BoolElt Default: true Precision : RngReSubElt Default: 0.00005 ``` Determines the threshold of the described ensemble of LDPC codes over the Gaussian channel, which is the critical value of the standard deviation above which there is a non-vanishing error probability (asymptotically). The ensemble can either be defined by two integers for $(v,c)$-regular LDPC codes, or by two density distributions $Sv$ and $Sc$, which are sequences of non-negative real numbers. The density distributions do not need to be normalized, though the first entry (corresponding to weight $1$ nodes in the Tanner graph) should always be zero. The computation proceeds by establishing lower and upper bounds on the threshold, then narrowing this range by repeatedly performing density evolution on the midpoint. If the threshold is approximately known then manually setting tight `Lower` and `Upper` bounds can reduce the length of the calculation. The speed with which these evolutions are computed depends on how fine the discretization is, controlled by the variable argument `Points`. If the threshold is needed to high levels of accuracy then an initial computation with fewer points is recommended to get a reduced searched range. The specific meaning of each variable argument is described below. `Lower` and `Upper` are real-valued bounds on the threshold, which (if tight) can help to reduce the search range and speed up the threshold determination. The validity of an input bound is verified before the search begins, and an error is returned if it is incorrect. `Points` and `MaxLLR` define the discretized basis of log likelihood ratios on which density evolution is performed, and have integer and real values resp. Specifically, the probability mass function is defined on the range $[-$`MaxLLR`$, \ldots,$`MaxLLR`$]$ on $2*Points+1$ discretized points. `MaxIterations` allows the user to set a finite limit of iterations that a density evolution should perform in determining the asymptotic behaviour at each channel parameter. Although this may help reduce the time of a computation, it should be kept in mind that the result may not be valid. `QuickCheck` defines the method by which the asymptotic behaviour at each channel parameter is identified. If set to `false`, then the probability density must evolve all the way to within an infinitesimal value of unity. When set to `true`, if the rate of change of the probability density is seen to be successively increasing then the asymptotic behaviour is assumed to go to unity. Empirically this method seems to give accurate results and so the default behaviour is `true`, however it has no theoretical justification. `Precision` is a real-valued parameter defining the precision to which the threshold should be determined. Setting the verbose mode `Code` prints out the bounds on the threshold as subsequent density evolutions narrow the search range. ### `DensityEvolutionGaussian(v, c, σ): RngIntElt, RngIntElt, FldReElt -> BoolElt` ### `DensityEvolutionGaussian(Sv, Sc, σ): [FldReElt], [FldReElt], FldReElt -> BoolElt` ```magma Points : RngIntElt Default: 500 MaxLLR : RngReSubElt Default: 25 MaxIterations: RngIntElt Default: Infinity() QuickCheck : BoolElt Default: true ``` Perform density evolution on the Gaussian channel using standard deviation $\sigma$ and determine the asymptotic behaviour for the given LDPC ensemble. The return value is boolean, where `true` indicates that $\sigma$ is below the threshold and the ensemble has error probability asymptotically tending to zero. See the description of `LDPCGaussianThreshold` for a description of the variable arguments. ### `GoodLDPCEnsemble(i): RngIntElt -> FldReElt, [FldReElt], [FldReElt]` Access a small database of density distributions defining good irregular LDPC ensembles. Returned is the published threshold of the ensemble over the Gaussian channel, along with the variable and check degree distributions. The input $i$ is a non-negative integer, which indexes the database (in no particular order). ### `Example: DE Gaussian (ex-4cd36e)` Since performing density evolution on a large number of discrete points is time consuming, it is normally better to first get an estimate with an easier computation. In this example a published value of the threshold of an ensemble (obtained using a different implementation) can be compared to the outputs from different levels of discretization. ```magma > thresh, Sv, Sc := GoodLDPCEnsemble(5); > R4 := RealField(4); > [R4| x : x in Sv]; [ 0.0000, 0.3001, 0.2839, 0.0000, 0.0000, 0.0000, 0.0000, 0.4159 ] > [R4| x : x in Sc]; [ 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.2292, 0.7708 ] > thresh; 0.949700000000000000000000000000 > time approx1 := LDPCGaussianThreshold(Sv, Sc: > Points := 500, Precision := 0.001); Time: 19.400 > approx1; 0.960976562499999964472863211995 > time approx2 := LDPCGaussianThreshold(Sv, Sc: > Points := 3000, > Lower := approx1-0.02, Upper := approx1+0.02); Time: 873.560 > approx2; 0.951210937499999964472863211995 ``` This estimate can now be used to narrow the search range of accurate density evolution. For the very long calculation the verbose mode is used to keep track of the progress of the calculation. ```magma > SetVerbose("Code", true); > time approx3 := LDPCGaussianThreshold(Sv, Sc: > Points := 5000, > Lower := approx2-0.005, Upper := approx2+0.0005); Determining the mapping matrix... ...mapping matrix obtained 19.10s Threshold Determination for LDPC code ensemble: c: (6):0.229190 (7):0.770810 v: (2):0.300130 (3):0.283950 (8):0.415920 will be found to precision 0.000050 Max LLR is 25.000000 distributed across 5000 points Asymptotic behaviour determination is: fast Beginning search with lb = 0.946211, ub = 0.951711 New Bounds: lb = 0.948961, ub = 0.951711 114.41s New Bounds: lb = 0.950336, ub = 0.951711 367.19s New Bounds: lb = 0.950336, ub = 0.951023 553.95s New Bounds: lb = 0.950336, ub = 0.950680 814.35s New Bounds: lb = 0.950508, ub = 0.950680 1261.46s New Bounds: lb = 0.950508, ub = 0.950594 1557.46s New Bounds: lb = 0.950508, ub = 0.950551 1891.52s Time: 2136.150 ``` The thresholds given in the database are published values taken from other implementations, and so are not guaranteed to match up exactly with the values obtained using Magma.