Low Density Lattice Codes

TL;DR

Low density lattice codes (LDLC) are a new class of lattice codes with sparse inverse generator matrices, enabling efficient decoding and near-capacity performance on the AWGN channel, demonstrated through large-scale simulations.

Contribution

This paper introduces LDLC, a novel lattice coding scheme with sparse inverse matrices, allowing linear-time decoding and near-capacity error performance.

Findings

Decoding scheme achieves near-capacity performance within ~0.5dB.

Efficient linear-time iterative decoding demonstrated.

Good error performance at block length of 100,000 symbols.

Abstract

Low density lattice codes (LDLC) are novel lattice codes that can be decoded efficiently and approach the capacity of the additive white Gaussian noise (AWGN) channel. In LDLC a codeword x is generated directly at the n-dimensional Euclidean space as a linear transformation of a corresponding integer message vector b, i.e., x = Gb, where H, the inverse of G, is restricted to be sparse. The fact that H is sparse is utilized to develop a linear-time iterative decoding scheme which attains, as demonstrated by simulations, good error performance within ~0.5dB from capacity at block length of n = 100,000 symbols. The paper also discusses convergence results and implementation considerations.

Figures

Equations

$$\mathbf{\boldsymbol{H}}=\left(\begin{array}[]{cccccc}0&-0.8&0&-0.5&1&0\\ 0.8&0&0&1&0&-0.5\\ 0&0.5&1&0&0.8&0\\ 0&0&-0.5&-0.8&0&1\\ 1&0&0&0&0.5&0.8\\ 0.5&-1&-0.8&0&0&0\end{array}\right)$$

$$f_{x_k|\underline{y}}(x_k|\underline{y})=C\cdot \sum _{\underline{l}\in \mathbf{G}\cap B}\delta (x_k-l_k)\cdot e^{-d^2(\underline{l},\underline{y})/2{\sigma }^2}$$

$${\tilde{p}}_j(x)=f_1\left(\frac{x}{h_1}\right)⊛\cdots f_{j-1}\left(\frac{x}{h_{j-1}}\right)⊛$$

$$⊛f_{j+1}\left(\frac{x}{h_{j+1}}\right)⊛\cdots \cdots ⊛f_r\left(\frac{x}{h_r}\right)$$

$$Q_j(x)=\sum _{i=-\infty }^{\infty }{p}_j\left(x-\frac{i}{h_j}\right)$$

$$\frac{1}{V_i^{(t+1)}}=\frac{1}{{\sigma }^2}+\sum _{\begin{array}{c}m=1\\ m\ne i\end{array}}^d\frac{h_m^2}{{\sum }_{\begin{array}{c}j=1\\ j\ne m\end{array}}^d{h}_j^2{V}_j^{(t)}}$$

$$\frac{1}{V_1^{(t+1)}}=\frac{h_2^2}{h_1^2{V}_1^{(t)}+h_3^2{V}_3^{(t)}}+\frac{h_3^2}{h_1^2{V}_1^{(t)}+h_2^2{V}_2^{(t)}}+\frac{1}{{\sigma }^2}$$

$$\frac{1}{V_2^{(t+1)}}=\frac{h_1^2}{h_2^2{V}_2^{(t)}+h_3^2{V}_3^{(t)}}+\frac{h_3^2}{h_1^2{V}_1^{(t)}+h_2^2{V}_2^{(t)}}+\frac{1}{{\sigma }^2}$$

$$\frac{1}{V_3^{(t+1)}}=\frac{h_1^2}{h_2^2{V}_2^{(t)}+h_3^2{V}_3^{(t)}}+\frac{h_2^2}{h_1^2{V}_1^{(t)}+h_3^2{V}_3^{(t)}}+\frac{1}{{\sigma }^2}$$

$$V_1^{(\infty )}\stackrel{\Delta }{=}\underset{t\to \infty }{\lim }{V}_1^{(t)}={\sigma }^2(1-\alpha ).$$

$$V_i^{(\infty )}\stackrel{\Delta }{=}\underset{t\to \infty }{\lim }{V}_i^{(t)}=0$$

$$0<\underset{t\to \infty }{\lim }\left|\frac{V_i^{(t)}-V_i^{(\infty )}}{{\alpha }^t}\right|<\infty$$

$$V_i^{(t)}\le {\sigma }^2{\alpha }^t$$

$${\underline{e}}^{(t+1)}\approx -\tilde{\mathbf{H}}\cdot {\underline{e}}^{(t)}$$

$${\underline{e}}^{(t+1)}\approx -\mathbf{F}\cdot {\underline{e}}^{(t)}+(1-\alpha )(\underline{y}-\mathbf{G}\underline{b})$$

$$\displaystyle F_{k,l}=\left\{\begin{array}[]{ll}\frac{H_{r,k}}{H_{r,l}}&\textrm{if }k\neq l\textrm{ and there exist a row }r\textrm{ of H}\\ &\textrm{for which }|H_{r,l}|=h_{1}\textrm{ and }H_{r,k}\neq 0\\ 0&\textrm{otherwise}\end{array}\right.$$

$${\underline{l}}^{(t+1)}=\mathbf{A}\cdot {\underline{l}}^{(t)}-{\underline{a}}^{(t)}-{\underline{c}}^{(t)}$$

$$a_{(k-1)d+i}^{(t)}=$$

$$=\frac{{\hat{V}}_{k,i}^{(t)}}{2}\left(\sum _{\begin{array}{c}l=1\\ l\ne i\end{array}}^d\sum _{\begin{array}{c}j=l+1\\ j\ne i\end{array}}^d\frac{{\left({\tilde{m}}_{k,l}^{(t)}-{\tilde{m}}_{k,j}^{(t)}\right)}^2}{{\tilde{V}}_{k,l}^{(t)}\cdot {\tilde{V}}_{k,j}^{(t)}}+\sum _{\begin{array}{c}l=1\\ l\ne i\end{array}}^d\frac{{\left({\tilde{m}}_{k,l}^{(t)}-y_k\right)}^2}{{\sigma }^2\cdot {\tilde{V}}_{k,l}^{(t)}}\right)$$

$$\mathbf{W}\stackrel{\Delta }{=}(d+1-\alpha )\mathbf{I}-2(1-\alpha )(\mathbf{I}+\mathbf{F}{)}^{-1}+$$

$$+(1-\alpha )(\mathbf{I}+\mathbf{F}{)}^{{-1}^T}\left((d-1{)}^2\mathbf{I}-{\mathbf{F}}^T\mathbf{F}\right)(\mathbf{I}+\mathbf{F}{)}^{-1}$$

$$s^{(t+1)}=(d-1{)}^2{s}^{(t)}-\sum _{i=1}^{nd}{a}_i^{(t)}$$

$${\tilde{s}}^{(t+1)}={\tilde{s}}^{(t)}-\frac{1}{(d-1{)}^{2t+2}}\sum _{i=1}^{nd}{a}_i^{(t)}$$

$${\tilde{s}}^{(t)}=-\sum _{j=0}^{t-1}\frac{{\sum }_{i=1}^{nd}{a}_i^{(j)}}{(d-1{)}^{2j+2}}$$

$$P^{(t)}=e^{s^{(t)}}=e^{(d-1{)}^{2t}\cdot {\tilde{s}}^{(t)}}$$

$$\underset{t\to \infty }{\lim }\frac{P_1^{(t)}}{P_0^{(t)}}=\frac{e^{-(d-1{)}^{2t}\cdot S_1}}{e^{-(d-1{)}^{2t}\cdot S_0}}=e^{(d-1{)}^{2t}\cdot (S_0-S_1)}$$

$$s^{(t+1)}=s^{(t)}-\sum _{i=1}^{2n}{a}_i^{(t)}$$

$$s^{(t)}=-\sum _{j=0}^{t-1}\sum _{i=1}^{2n}{a}_i^{(j)}$$

$$\underset{t\to \infty }{\lim }\frac{P_1^{(t)}}{P_0^{(t)}}=\frac{e^{-S_{a1}\cdot t}}{e^{-S_{a0}\cdot t}}=e^{(S_{a0}-S_{a1})\cdot t}$$

$$S_a\stackrel{\Delta }{=}\underset{j\to \infty }{\lim }\sum _{i=1}^{2n}{a}_i^{(j)}=\frac{1}{2{\sigma }^2}(\mathbf{G}\underline{b}-\underline{y}{)}^T\mathbf{W}(\mathbf{G}\underline{b}-\underline{y})$$

$${\underline{x}}^{(t)}=\underline{\tilde{b}}-\tilde{\mathbf{H}}\cdot {\underline{x}}^{(t-1)}$$

Full text

I Introduction

If we take a look at the evolution of codes for binary or finite alphabet channels, it was first shown [1] that channel capacity can be achieved with long random codewords. Then, it was found out [2] that capacity can be achieved via a simpler structure of linear codes. Then, specific families of linear codes were found that are practical and have good minimum Hamming distance (e.g. convolutional codes, cyclic block codes, specific cyclic codes such as BCH and Reed-Solomon codes [4]). Later, capacity achieving schemes were found, which have special structures that allow efficient iterative decoding, such as low-density parity-check (LDPC) codes [5] or turbo codes [6].

If we now take a similar look at continuous alphabet codes for the additive white Gaussian noise (AWGN) channel, it was first shown [3] that codes with long random Gaussian codewords can achieve capacity. Later, it was shown that lattice codes can also achieve capacity ([7] – [12]). Lattice codes are clearly the Euclidean space analogue of linear codes. Similarly to binary codes, we could expect that specific practical lattice codes will then be developed. However, there was almost no further progress from that point. Specific lattice codes that were found were based on fixed dimensional classical lattices [19] or based on algebraic error correcting codes [13][14], but no significant effort was made in designing lattice codes directly in the Euclidean space or in finding specific capacity achieving lattice codes. Practical coding schemes for the AWGN channel were based on finite alphabet codes.

In [15], “signal codes” were introduced. These are lattice codes, designed directly in the Euclidean space, where the information sequence of integers $i_{n}$, $n=1,2,...$ is encoded by convolving it with a fixed signal pattern $g_{n}$, $n=1,2,...d$. Signal codes are clearly analogous to convolutional codes, and in particular can work at the AWGN channel cutoff rate with simple sequential decoders. In [16] it is also demonstrated that signal codes can work near the AWGN channel capacity with more elaborated bi-directional decoders. Thus, signal codes provided the first step toward finding effective lattice codes with practical decoders.

Inspired by LDPC codes and in the quest of finding practical capacity achieving lattice codes, we propose in this work “Low Density Lattice Codes” (LDLC). We show that these codes can approach the AWGN channel capacity with iterative decoders whose complexity is linear in block length. In recent years several schemes were proposed for using LDPC over continuous valued channels by either multilevel coding [18] or by non-binary alphabet (e.g. [17]). Unlike these LDPC based schemes, in LDLC both the encoder and the channel use the same real algebra which is natural for the continuous-valued AWGN channel. This feature also simplifies the convergence analysis of the iterative decoder.

The outline of this paper is as follows. Low density lattice codes are first defined in Section II. The iterative decoder is then presented in Section III, followed by convergence analysis of the decoder in Section IV. Then, Section V describes how to choose the LDLC code parameters, and Section VI discusses implementation considerations. The computational complexity of the decoder is then discussed in Section VII, followed by a brief description of encoding and shaping in Section VIII. Simulation results are finally presented in Section IX.

II Basic Definitions and Properties

II-A Lattices and Lattice Codes

An $n$ dimensional lattice in $\mathbb{R}^{m}$ is defined as the set of all linear combinations of a given basis of $n$ linearly independent vectors in $\mathbb{R}^{m}$ with integer coefficients. The matrix $\boldsymbol{G}$, whose columns are the basis vectors, is called a generator matrix of the lattice. Every lattice point is therefore of the form $\underline{x}=\boldsymbol{G}\underline{b}$, where $\underline{b}$ is an $n$-dimensional vector of integers. The Voronoi cell of a lattice point is defined as the set of all points that are closer to this point than to any other lattice point. The Voronoi cells of all lattice points are congruent, and for square $\boldsymbol{G}$ the volume of the Voronoi cell is equal to $det(\boldsymbol{G})$. In the sequel $\boldsymbol{G}$ will be used to denote both the lattice and its generator matrix.

A lattice code of dimension $n$ is defined by a (possibly shifted) lattice $\boldsymbol{G}$ in $\mathbb{R}^{m}$ and a shaping region $B\subset\mathbb{R}^{m}$, where the codewords are all the lattice points that lie within the shaping region $B$. Denote the number of these codewords by $N$. The average transmitted power (per channel use, or per symbol) is the average energy of all codewords, divided by the codeword length $m$. The information rate (in bits/symbol) is $log_{2}(N)/m$.

When using a lattice code for the AWGN channel with power limit $P$ and noise variance $\sigma^{2}$, the maximal information rate is limited by the capacity $\frac{1}{2}\log_{2}(1+\frac{P}{\sigma^{2}})$. Poltyrev [20] considered the AWGN channel without restrictions. If there is no power restriction, code rate is a meaningless measure, since it can be increased without limit. Instead, it was suggested in [20] to use the measure of constellation density, leading to a generalized definition of the capacity as the maximal possible codeword density that can be recovered reliably. When applied to lattices, the generalized capacity implies that there exists a lattice $\boldsymbol{G}$ of high enough dimension $n$ that enables transmission with arbitrary small error probability, if and only if $\sigma^{2}<\frac{\sqrt[n]{|det(\boldsymbol{G})|^{2}}}{2\pi e}$. A lattice that achieves the generalized capacity of the AWGN channel without restrictions, also achieves the channel capacity of the power constrained AWGN channel, with a properly chosen spherical shaping region (see also [12]).

In the rest of this work we shall concentrate on the lattice design and the lattice decoding algorithm, and not on the shaping region or shaping algorithms. We shall use lattices with $det(\boldsymbol{G})=1$, where analysis and simulations will be carried for the AWGN channel without restrictions. A capacity achieving lattice will have small error probability for noise variance $\sigma^{2}$ which is close to the theoretical limit $\frac{1}{2\pi e}$.

II-B Syndrome and Parity Check Matrix for Lattice Codes

A binary $(n,k)$ error correcting code is defined by its $n\times k$ binary generator matrix $\boldsymbol{G}$. A binary information vector $\underline{b}$ with dimension $k$ is encoded by $\underline{x}=\boldsymbol{G}\underline{b}$, where calculations are performed in the finite field GF(2). The parity check matrix $\boldsymbol{H}$ is an $(n-k)\times n$ matrix such that $\underline{x}$ is a codeword if and only if $\boldsymbol{H}\underline{x}=\underline{0}$. The input to the decoder is the noisy codeword $\underline{y}=\underline{x}+\underline{e}$, where $\underline{e}$ is the error sequence and addition is done in the finite field. The decoder typically starts by calculating the syndrome $\underline{s}=\boldsymbol{H}\underline{y}=\boldsymbol{H}(\underline{x}+\underline{e})=\boldsymbol{H}\underline{e}$ which depends only on the noise sequence and not on the transmitted codeword.

We would now like to extend the definitions of the parity check matrix and the syndrome to lattice codes. An $n$-dimensional lattice code is defined by its $n\times n$ lattice generator matrix $\boldsymbol{G}$ (throughout this paper we assume that $\boldsymbol{G}$ is square, but the results are easily extended to the non-square case). Every codeword is of the form $\underline{x}=\boldsymbol{G}\underline{b}$, where $\underline{b}$ is a vector of integers. Therefore, $\boldsymbol{G}^{-1}\underline{x}$ is a vector of integers for every codeword $\underline{x}$. We define the parity check matrix for the lattice code as $\boldsymbol{H}\stackrel{{\scriptstyle\Delta}}{{=}}\boldsymbol{G}^{-1}$. Given a noisy codeword $\underline{y}=\underline{x}+\underline{w}$ (where $w$ is the additive noise vector, e.g. AWGN, added by real arithmetic), we can then define the syndrome as $\underline{s}\stackrel{{\scriptstyle\Delta}}{{=}}frac\{\boldsymbol{H}\underline{y}\}$, where $frac\{x\}$ is the fractional part of $x$, defined as $frac\{x\}=x-\left\lfloor x\right\rceil$, where $\left\lfloor x\right\rceil$ denotes the nearest integer to $x$. The syndrome $\underline{s}$ will be zero if and only if $\underline{y}$ is a lattice point, since $\boldsymbol{H}\underline{y}$ will then be a vector of integers with zero fractional part. For a noisy codeword, the syndrome will equal $\underline{s}=frac\{\boldsymbol{H}\underline{y}\}=frac\{\boldsymbol{H}(\underline{x}+\underline{w})\}=frac\{\boldsymbol{H}\underline{w}\}$ and therefore will depend only on the noise sequence and not on the transmitted codeword, as desired.

Note that the above definitions of the syndrome and parity check matrix for lattice codes are consistent with the definitions of the dual lattice and the dual code[19]: the dual lattice of a lattice $\boldsymbol{G}$ is defined as the lattice with generator matrix $\boldsymbol{H}=\boldsymbol{G}^{-1}$, where for binary codes, the dual code of $\boldsymbol{G}$ is defined as the code whose generator matrix is $\boldsymbol{H}$, the parity check matrix of $\boldsymbol{G}$.

II-C Low Density Lattice Codes

We shall now turn to the definition of the codes proposed in this paper - low density lattice codes (LDLC).

Definition 1** (LDLC)**

An $n$ dimensional LDLC is an $n$-dimensional lattice code with a non-singular lattice generator matrix $\boldsymbol{G}$ satisfying $|det(\boldsymbol{G})|=1$, for which the parity check matrix $\boldsymbol{H}=\boldsymbol{G}^{-1}$ is sparse. The $i$’th row degree $r_{i}$, $i=1,2,...n$ is defined as the number of nonzero elements in row $i$ of $\boldsymbol{H}$, and the $i$’th column degree $c_{i}$, $i=1,2,...n$ is defined as the number of nonzero elements in column $i$ of $\boldsymbol{H}$.

Note that in binary LDPC codes, the code is completely defined by the locations of the nonzero elements of $\boldsymbol{H}$. In LDLC there is another degree of freedom since we also have to choose the values of the nonzero elements of $\boldsymbol{H}$.

Definition 2** (regular LDLC)**

An $n$ dimensional LDLC is regular if all the row degrees and column degrees of the parity check matrix are equal to a common degree $d$.

Definition 3** (magic square LDLC)**

An $n$ dimensional regular LDLC with degree $d$ is called “magic square LDLC” if every row and column of the parity check matrix $\boldsymbol{H}$ has the same $d$ nonzero values, except for a possible change of order and random signs. The sorted sequence of these $d$ values $h_{1}\geq h_{2}\geq...\geq h_{d}>0$ will be referred to as the generating sequence of the magic square LDLC.

For example, the matrix

[TABLE]

is a parity check matrix of a magic square LDLC with lattice dimension $n=6$, degree $d=3$ and generating sequence $\{1,0.8,0.5\}$. This $\boldsymbol{H}$ should be further normalized by the constant $\sqrt[n]{|det(\boldsymbol{H})|}$ in order to have $|det(\boldsymbol{H})|=|det(\boldsymbol{G})|=1$, as required by Definition 1.

The bipartite graph of an LDLC is defined similarly to LDPC codes: it is a graph with variable nodes at one side and check nodes at the other side. Each variable node corresponds to a single element of the codeword $\underline{x}=\boldsymbol{G}\underline{b}$. Each check node corresponds to a check equation (a row of $\boldsymbol{H}$). A check equation is of the form $\sum_{k}h_{k}x_{i_{k}}=integer$, where $i_{k}$ denotes the locations of the nonzero elements at the appropriate row of $\boldsymbol{H}$, $h_{k}$ are the values of $\boldsymbol{H}$ at these locations and the integer at the right hand side is unknown. An edge connects check node $i$ and variable node $j$ if and only if $H_{i,j}\neq 0$. This edge is assigned the value $H_{i,j}$. Figure 1 illustrates the bi-partite graph of a magic square LDLC with degree 3. In the figure, every variable node $x_{k}$ is also associated with its noisy channel observation $y_{k}$.

Finally, a $k$-loop is defined as a loop in the bipartite graph that consists of $k$ edges. A bipartite graph, in general, can only contain loops with even length. Also, a 2-loop, which consists of two parallel edges that originate from the same variable node to the same check node, is not possible by the definition of the graph. However, longer loops are certainly possible. For example, a 4-loop exists when two variable nodes are both connected to the same pair of check nodes.

III Iterative Decoding for the AWGN Channel

Assume that the codeword $\underline{x}=\boldsymbol{G}\underline{b}$ was transmitted, where $\underline{b}$ is a vector of integers. We observe the noisy codeword $\underline{y}=\underline{x}+\underline{w}$, where $\underline{w}$ is a vector of i.i.d Gaussian noise samples with common variance $\sigma^{2}$, and we need to estimate the integer valued vector $\underline{b}$. The maximum likelihood (ML) estimator is then $\hat{\underline{b}}=\arg\mathop{\min}\limits_{\underline{b}}||\underline{y}-\boldsymbol{G}\underline{b}||^{2}$.

Our decoder will not estimate directly the integer vector $\underline{b}$. Instead, it will estimate the probability density function (PDF) of the codeword vector $\underline{x}$. Furthermore, instead of calculating the $n$-dimensional PDF of the whole vector $\underline{x}$, we shall calculate the $n$ one-dimensional PDF’s for each of the components $x_{k}$ of this vector (conditioned on the whole observation vector $\underline{y}$). In appendix A it is shown that $f_{x_{k}|\underline{y}}(x_{k}|\underline{y})$ is a weighted sum of Dirac delta functions:

[TABLE]

where $\underline{l}$ is a lattice point (vector), $l_{k}$ is its $k$-th component, $C$ is a constant independent of $x_{k}$ and $d(\underline{l},\underline{y})$ is the Euclidean distance between $\underline{l}$ and $\underline{y}$. Direct evaluation of (1) is not practical, so our decoder will try to estimate $f_{x_{k}|\underline{y}}(x_{k}|\underline{y})$ (or at least approximate it) in an iterative manner.

Our decoder will decode to the infinite lattice, thus ignoring the shaping region boundaries. This approximate decoding method is no longer exact maximum likelihood decoding, and is usually denoted “lattice decoding” [12].

The calculation of $f_{x_{k}|\underline{y}}(x_{k}|\underline{y})$ is involved since the components $x_{k}$ are not independent random variables (RV’s), because $\underline{x}$ is restricted to be a lattice point. Following [5] we use a “trick” - we assume that the $x_{k}$’s are independent, but add a condition that assures that $\underline{x}$ is a lattice point. Specifically, define $\underline{s}\stackrel{{\scriptstyle\Delta}}{{=}}\boldsymbol{H}\cdot\underline{x}$. Restricting $\underline{x}$ to be a lattice point is equivalent to restricting $\underline{s}\in\mathbb{Z}^{n}$. Therefore, instead of calculating $f_{x_{k}|\underline{y}}(x_{k}|\underline{y})$ under the assumption that $\underline{x}$ is a lattice point, we can calculate $f_{x_{k}|\underline{y}}(x_{k}|\underline{y},\underline{s}\in\mathbb{Z}^{n})$ and assume that the $x_{k}$ are independent and identically distributed (i.i.d) with a continuous PDF (that does not include Dirac delta functions). It still remains to set $f_{x_{k}}(x_{k})$, the PDF of $x_{k}$. Under the i.i.d assumption, the PDF of the codeword $\underline{x}$ is $f_{\underline{x}}(\underline{x})=\prod_{k=1}^{n}f_{x_{k}}(x_{k})$. As shown in Appendix B, the value of $f_{\underline{x}}(\underline{x})$ is not important at values of $\underline{x}$ which are not lattice points, but at a lattice point it should be proportional to the probability of using this lattice point. Since we assume that all lattice points are used equally likely, $f_{\underline{x}}(\underline{x})$ must have the same value at all lattice points. A reasonable choice for $f_{x_{k}}(x_{k})$ is then to use a uniform distribution such that $\underline{x}$ will be uniformly distributed in an $n$-dimensional cube. For an exact ML decoder (that takes into account the boundaries of the shaping region), it is enough to choose the range of $f_{x_{k}}(x_{k})$ such that this cube will contain the shaping region. For our decoder, that performs lattice decoding, we should set the range of $f_{x_{k}}(x_{k})$ large enough such that the resulting cube will include all the lattice points which are likely to be decoded. The derivation of the iterative decoder shows that this range can be set as large as needed without affecting the complexity of the decoder.

The derivation in [5] further imposed the tree assumption. In order to understand the tree assumption, it is useful to define the tier diagram, which is shown in Figure 2 for a regular LDLC with degree 3. Each vertical line corresponds to a check equation. The tier 1 nodes of $x_{1}$ are all the elements $x_{k}$ that take place in a check equation with $x_{1}$. The tier 2 nodes of $x_{1}$ are all the elements that take place in check equations with the tier 1 elements of $x_{1}$, and so on. The tree assumption assumes that all the tree elements are distinct (i.e. no element appears in different tiers or twice in the same tier). This assumption simplifies the derivation, but in general, does not hold in practice, so our iterative algorithm is not guaranteed to converge to the exact solution (1) (see Section IV).

The detailed derivation of the iterative decoder (using the above “trick” and the tree assumption) is given in Appendix C. In Section III-A below we present the final resulting algorithm. This iterative algorithm can also be explained by intuitive arguments, described after the algorithm specification.

III-A The Iterative Decoding Algorithm

The iterative algorithm is most conveniently represented by using a message passing scheme over the bipartite graph of the code, similarly to LDPC codes. The basic difference is that in LDPC codes the messages are scalar values (e.g. the log likelihood ratio of a bit), where for LDLC the messages are real functions over the interval $(-\infty,\infty)$. As in LDPC, in each iteration the check nodes send messages to the variable nodes along the edges of the bipartite graph and vice versa. The messages sent by the check nodes are periodic extensions of PDF’s. The messages sent by the variable nodes are PDF’s.

LDLC iterative decoding algorithm:

Denote the variable nodes by $x_{1},x_{2},...,x_{n}$ and the check nodes by $c_{1},c_{2},...c_{n}$.

• •

Initialization: each variable node $x_{k}$ sends to all its check nodes the message $f_{k}^{(0)}(x)=\frac{1}{\sqrt{2\pi\sigma^{2}}}e^{-\frac{(y_{k}-x)^{2}}{2\sigma^{2}}}$.

• •

Basic iteration - check node message: Each check node sends a (different) message to each of the variable nodes that are connected to it. For a specific check node denote (without loss of generality) the appropriate check equation by $\sum_{l=1}^{r}h_{l}x_{m_{l}}=integer$, where $x_{m_{l}}$, $l=1,2...r$ are the variable nodes that are connected to this check node (and $r$ is the appropriate row degree of $\boldsymbol{H}$). Denote by $f_{l}(x)$, $l=1,2...r$, the message that was sent to this check node by variable node $x_{m_{l}}$ in the previous half-iteration. The message that the check node transmits back to variable node $x_{m_{j}}$ is calculated in three basic steps.

The convolution step - all messages, except $f_{j}(x)$, are convolved (after expanding each $f_{l}(x)$ by $h_{l}$):

[TABLE] 2. 2.

The stretching step - The result is stretched by $(-h_{j})$ to $p_{j}(x)=\tilde{p}_{j}(-h_{j}x)$ 3. 3.

The periodic extension step - The result is extended to a periodic function with period $1/|h_{j}|$:

[TABLE]

The function $Q_{j}(x)$ is the message that is finally sent to variable node $x_{m_{j}}$.

• •

Basic iteration - variable node message: Each variable node sends a (different) message to each of the check nodes that are connected to it. For a specific variable node $x_{k}$, assume that it is connected to check nodes $c_{m_{1}},c_{m_{2}},...c_{m_{e}}$, where $e$ is the appropriate column degree of $\boldsymbol{H}$. Denote by $Q_{l}(x)$, $l=1,2,...e$, the message that was sent from check node $c_{m_{l}}$ to this variable node in the previous half-iteration. The message that is sent back to check node $c_{m_{j}}$ is calculated in two basic steps:

The product step: $\tilde{f}_{j}(x)=e^{-\frac{(y_{k}-x)^{2}}{2\sigma^{2}}}\prod_{\begin{subarray}{c}l=1\\ l\neq j\end{subarray}}^{e}Q_{l}(x)$ 2. 2.

The normalization step: $f_{j}(x)=\frac{\tilde{f}_{j}(x)}{\int_{-\infty}^{\infty}\tilde{f}_{j}(x)dx}$

This basic iteration is repeated for the desired number of iterations.

• •

Final decision: After finishing the iterations, we want to estimate the integer information vector $\underline{b}$. First, we estimate the final PDF’s of the codeword elements $x_{k}$, $k=1,2,...n$, by calculating the variable node messages at the last iteration without omitting any check node message in the product step: $\tilde{f}^{(k)}_{final}(x)=e^{-\frac{(y_{k}-x)^{2}}{2\sigma^{2}}}\prod_{l=1}^{e}Q_{l}(x)$. Then, we estimate each $x_{k}$ by finding the peak of its PDF: $\hat{x_{k}}=arg\max_{x}\tilde{f}^{(k)}_{final}(x)$. Finally, we estimate $\underline{b}$ as $\underline{\hat{b}}=\left\lfloor\boldsymbol{H}\underline{\hat{x}}\right\rceil$.

The operation of the iterative algorithm can be intuitively explained as follows. The check node operation is equivalent to calculating the PDF of $x_{m_{j}}$ from the PDF’s of $x_{m_{i}}$, $i=1,2,...,j-1,j+1,...r$, given that $\sum_{l=1}^{r}h_{l}x_{m_{l}}=integer$, and assuming that $x_{m_{i}}$ are independent. Extracting $x_{m_{j}}$ from the check equation, we get $x_{m_{j}}=\frac{1}{h_{j}}(integer-\sum_{\begin{subarray}{c}l=1\\ l\neq j\end{subarray}}^{r}h_{l}x_{m_{l}})$. Since the PDF of a sum of independent random variables is the convolution of the corresponding PDF’s, equation (1) and the stretching step that follows it simply calculate the PDF of $x_{m_{j}}$, assuming that the integer at the right hand side of the check equation is zero. The result is then periodically extended such that a properly shifted copy exists for every possible value of this (unknown) integer. The variable node gets such a message from all the check equations that involve the corresponding variable. The check node messages and the channel PDF are treated as independent sources of information on the variable, so they are multiplied all together.

Note that the periodic extension step at the check nodes is equivalent to a convolution with an infinite impulse train. With this observation, the operation of the variable nodes is completely analogous to that of the check nodes: the variable nodes multiply the incoming messages by the channel PDF, where the check nodes convolve the incoming messages with an impulse train, which can be regarded as a generalized “integer PDF”.

In the above formulation, the integer information vector $\underline{b}$ is recovered from the PDF’s of the codeword elements $x_{k}$. An alternative approach is to calculate the PDF of each integer element $b_{m}$ directly as the PDF of the left hand side of the appropriate check equation. Using the tree assumption, this can be done by simply calculating the convolution $\tilde{p}(x)$ as in (1), but this time without omitting any PDF, i.e. all the received variable node messages are convolved. Then, the integer $b_{m}$ is determined by $\hat{b}_{m}=arg\max_{j\in\mathbb{Z}}\tilde{p}(j)$.

Figure 3 shows an example for a regular LDLC with degree $d=5$. The figure shows all the signals that are involved in generating a variable node message for a certain variable node. The top signal is the channel Gaussian, centered around the noisy observation of the variable. The next 4 signals are the periodically extended PDF’s that arrived from the check nodes, and the bottom signal is the product of all the 5 signals. It can be seen that each periodic signal has a different period, according to the relevant coefficient of $\boldsymbol{H}$. Also, the signals with larger period have larger variance. This diversity resolves all the ambiguities such that the multiplication result (bottom plot) remains with a single peak. We expect the iterative algorithm to converge to a solution where a single peak will remain at each PDF, located at the desired value and narrow enough to estimate the information.

IV Convergence

IV-A The Gaussian Mixture Model

Interestingly, for LDLC we can come up with a convergence analysis that in many respects is more specific than the similar analysis for LDPC.

We start by introducing basic claims about Gaussian PDF’s. Denote $G_{m,V}(x)=\frac{1}{\sqrt{2\pi V}}e^{-\frac{(x-m)^{2}}{2V}}$.

Claim 1** (convolution of Gaussians)**

The convolution of $n$ Gaussians with mean values $m_{1},m_{2},...,m_{n}$ and variances $V_{1},V_{2},...,V_{n}$, respectively, is a Gaussian with mean $m_{1}+m_{2}+...+m_{n}$ and variance $V_{1}+V_{2}+...+V_{n}$.

Proof:

See [21]. ∎

Claim 2** (product of $n$ Gaussians)**

Let $G_{m_{1},V_{1}}(x)$, $G_{m_{2},V_{2}}(x)$,…,$G_{m_{n},V_{n}}(x)$ be $n$ Gaussians with mean values $m_{1},m_{2},...,m_{n}$ and variances $V_{1},V_{2},...,V_{n}$ respectively. Then, the product of these Gaussians is a scaled Gaussian: $\prod_{i=1}^{n}G_{m_{i},V_{i}}(x)=\hat{A}\cdot G_{\hat{m},\hat{V}}(x)$, where $\frac{1}{\hat{V}}=\sum_{i=1}^{n}\frac{1}{V_{i}}$, $\hat{m}=\frac{\sum_{i=1}^{n}m_{i}V_{i}^{-1}}{\sum_{i=1}^{n}V_{i}^{-1}}$, and $\hat{A}=\frac{1}{\sqrt{(2\pi)^{n-1}\hat{V}^{-1}\prod_{k=1}^{n}V_{k}}}\cdot e^{-\frac{\hat{V}}{2}\sum_{i=1}^{n}\sum_{j=i+1}^{n}\frac{(m_{i}-m_{j})^{2}}{V_{i}\cdot V_{j}}}$.

Proof:

By straightforward mathematical manipulations. ∎

The reason that we are interested in the properties of Gaussian PDF’s lies in the following lemma.

Lemma 1

Each message that is exchanged between the check nodes and variable nodes in the LDLC decoding algorithm (i.e. $Q_{j}(x)$ and $f_{j}(x)$), at every iteration, can be expressed as a Gaussian mixture of the form $M(x)=\sum_{j=1}^{\infty}A_{j}G_{m_{j},V_{j}}(x)$.

Proof:

By induction: The initial messages are Gaussians, and the basic operations of the iterative decoder preserve the Gaussian mixture nature of Gaussian mixture inputs (convolution and multiplication preserve the Gaussian nature according to claims 1 and 2, stretching, expanding and shifting preserve it by the definition of a Gaussian, and periodic extension transforms a single Gaussian to a mixture and a mixture to a mixture). ∎

Convergence analysis should therefore analyze the convergence of the variances, mean values and amplitudes of the Gaussians in each mixture.

IV-B Convergence of the Variances

We shall now analyze the behavior of the variances, and start with the following lemma.

Lemma 2

For both variable node messages and check node messages, all the Gaussians that take place in the same mixture have the same variance.

Proof:

By induction. The initial variable node messages are single element mixtures so the claim obviously holds. Assume now that all the variable node messages at iteration $t$ are mixtures where all the Gaussians that take place in the same mixture have the same variance. In the convolution step (1), each variable node message is first expanded. All Gaussians in the expanded mixture will still have the same variance, since the whole mixture is expanded together. Then, $d-1$ expanded Gaussian mixtures are convolved. In the resulting mixture, each Gaussian will be the result of convolving $d-1$ single Gaussians, one from each mixture. According to claim 1, all the Gaussians in the convolution result will have the same variance, which will equal the sum of the $d-1$ variances of the expanded messages. The stretching and periodic extension (3) do not change the equal variance property, so it holds for the final check node messages. The variable nodes multiply $d-1$ check node messages. Each Gaussian in the resulting mixture is a product of $d-1$ single Gaussians, one from each mixture, and the channel noise Gaussian. According to claim 2, they will all have the same variance. The final normalization step does not change the variances so the equal variance property is kept for the final variable node messages at iteration $t+1$. ∎

Until this point we did not impose any restrictions on the LDLC. From now on, we shall restrict ourselves to magic square regular LDLC (see Definition 3). The basic iterative equations that relate the variances at iteration $t+1$ to the variances at iteration $t$ are summarized in the following two lemmas.

Lemma 3

For magic square LDLC, variable node messages that are sent at the same iteration along edges with the same absolute value have the same variance.

Proof:

See Appendix D. ∎

Lemma 4

For magic square LDLC with degree $d$, denote the variance of the messages that are sent at iteration $t$ along edges with weight $\pm h_{l}$ by $V^{(t)}_{l}$. The variance values $V^{(t)}_{1},V^{(t)}_{2},...,V^{(t)}_{d}$ obey the following recursion:

[TABLE]

for $i=1,2,...d$, with initial conditions $V^{(0)}_{1}=V^{(0)}_{2}=...=V^{(0)}_{d}=\sigma^{2}$.

Proof:

See Appendix D. ∎

For illustration, the recursion for the case $d=3$ is:

[TABLE]

[TABLE]

[TABLE]

The lemmas above are used to prove the following theorem regarding the convergence of the variances.

Theorem 1

For a magic square LDLC with degree $d$ and generating sequence $h_{1}\geq h_{2}\geq...\geq h_{d}>0$, define $\alpha\stackrel{{\scriptstyle\Delta}}{{=}}\frac{\sum_{i=2}^{d}h_{i}^{2}}{h_{1}^{2}}$. Assume that $\alpha<1$. Then:

The first variance approaches a constant value of $\sigma^{2}(1-\alpha)$, where $\sigma^{2}$ is the channel noise variance:

[TABLE] 2. 2.

The other variances approach zero:

[TABLE]

for $i=2,3..d$. 3. 3.

The asymptotic convergence rate of all variances is exponential:

[TABLE]

for $i=1,2..d$. 4. 4.

The zero approaching variances are upper bounded by the decaying exponential $\sigma^{2}\alpha^{t}$:

[TABLE]

for $i=2,3..d$ and $t\geq 0$.

Proof:

See Appendix D. ∎

If $\alpha\geq 1$, the variances may still converge, but convergence rate may be as slow as $o(1/t)$, as illustrated in Appendix D.

Convergence of the variances to zero implies that the Gaussians approach impulses. This is a desired property of the decoder, since the exact PDF that we want to calculate is indeed a weighted sum of impulses (see (1)). It can be seen that by designing a code with $\alpha<1$, i.e. $h_{1}^{2}>\sum_{i=2}^{d}h_{i}^{2}$, one variance approaches a constant (and not zero). However, all the other variances approach zero, where all variances converge in an exponential rate. This will be the preferred mode because the information can be recovered even if a single variance does not decay to zero, where exponential convergence is certainly preferred over slow $1/t$ convergence. Therefore, from now on we shall restrict our analysis to magic square LDLC with $\alpha<1$.

Theorem 1 shows that every iteration, each variable node will generate $d-1$ messages with variances that approach zero, and a single message with variance that approaches a constant. The message with nonzero variance will be transmitted along the edge with largest weight (i.e. $h_{1}$). However, from the derivation of Appendix D it can be seen that the opposite happens for the check nodes: each check node will generate $d-1$ messages with variances that approach a constant, and a single message with variance that approaches zero. The check node message with zero approaching variance will be transmitted along the edge with largest weight.

IV-C Convergence of the Mean Values

The reason that the messages are mixtures and not single Gaussians lies in the periodic extension step (3) at the check nodes, and every Gaussian at the output of this step can be related to a single index of the infinite sum. Therefore, we can label each Gaussian at iteration $t$ with a list of all the indices that were used in (3) during its creation process in iterations $1,2,...t$.

Definition 4** (label of a Gaussian)**

The label of a Gaussian consists of a sequence of triplets of the form $\{t,c,i\}$, where $t$ is an iteration index, $c$ is a check node index and $i$ is an integer. The labels are initialized to the empty sequence. Then, the labels are updated along each iteration according to the following update rules:

In the periodic extension step (3), each Gaussian in the output periodic mixture is assigned the label of the specific Gaussian of $p_{j}(x)$ that generated it, concatenated with a single triplet $\{t,c,i\}$, where $t$ is the current iteration index, $c$ is the check node index and $i$ is the index in the infinite sum of (3) that corresponds to this Gaussian. 2. 2.

In the convolution step and the product step, each Gaussian in the output mixture is assigned a label that equals the concatenation of all the labels of the specific Gaussians in the input messages that formed this Gaussian. 3. 3.

The stretching and normalization steps do not alter the label of each Gaussian: Each Gaussian in the stretched/normalized mixture inherits the label of the appropriate Gaussian in the original mixture.

Definition 5** (a consistent Gaussian)**

A Gaussian in a mixture is called “[$t_{a}$, $t_{b}$] consistent” if its label contains no contradictions for iterations $t_{a}$ to $t_{b}$, i.e. for every pair of triplets $\{t_{1},c_{1},i_{1}\}$, $\{t_{2},c_{2},i_{2}\}$ such that $t_{a}\leq t_{1},t_{2}\leq t_{b}$, if $c_{1}=c_{2}$ then $i_{1}=i_{2}$. A [[math], $\infty$] consistent Gaussian will be simply called a consistent Gaussian.

We can relate every consistent Gaussian to a unique integer vector $\underline{b}\in\mathbb{Z}^{n}$, which holds the $n$ integers used in the $n$ check nodes. Since in the periodic extension step (3) the sum is taken over all integers, a consistent Gaussian exists in each variable node message for every possible integer valued vector $\underline{b}\in\mathbb{Z}^{n}$. We shall see later that this consistent Gaussian corresponds to the lattice point $\boldsymbol{G}\underline{b}$.

According to Theorem 1, if we choose the nonzero values of $\boldsymbol{H}$ such that $\alpha<1$, every variable node generates $d-1$ messages with variances approaching zero and a single message with variance that approaches a constant. We shall refer to these messages as “narrow” messages and “wide” messages, respectively. For a given integer valued vector $\underline{b}$, we shall concentrate on the consistent Gaussians that relate to $\underline{b}$ in all the $nd$ variable node messages that are generated in each iteration (a single Gaussian in each message). The following lemmas summarize the asymptotic behavior of the mean values of these consistent Gaussians for the narrow messages.

Lemma 5

For a magic square LDLC with degree $d$ and $\alpha<1$, consider the $d-1$ narrow messages that are sent from a specific variable node. Consider further a single Gaussian in each message, which is the consistent Gaussian that relates to a given integer vector $\underline{b}$. Asymptotically, the mean values of these $d-1$ Gaussians become equal.

Proof:

See Appendix E. ∎

Lemma 6

For a magic square LDLC with dimension $n$, degree $d$ and $\alpha<1$, consider only consistent Gaussians that relate to a given integer vector $\underline{b}$ and belong to narrow messages. Denote the common mean value of the $d-1$ such Gaussians that are sent from variable node $i$ at iteration $t$ by $m^{(t)}_{i}$, and arrange all these mean values in a column vector $\underline{m}^{(t)}$ of dimension $n$. Define the error vector $\underline{e}^{(t)}\stackrel{{\scriptstyle\Delta}}{{=}}\underline{m}^{(t)}-\underline{x}$, where $\underline{x}=\boldsymbol{G}\underline{b}$ is the lattice point that corresponds to $\underline{b}$. Then, for large $t$, $\underline{e}^{(t)}$ satisfies:

[TABLE]

where $\tilde{\boldsymbol{H}}$ is derived from $\boldsymbol{H}$ by permuting the rows such that the $\pm h_{1}$ elements will be placed on the diagonal, dividing each row by the appropriate diagonal element ($h_{1}$ or $-h_{1}$), and then nullifying the diagonal.

Proof:

See Appendix E. ∎

We can now state the following theorem, which describes the conditions for convergence and the steady state value of the mean values of the consistent Gaussians of the narrow variable node messages.

Theorem 2

For a magic square LDLC with $\alpha<1$, the mean values of the consistent Gaussians of the narrow variable node messages that relate to a given integer vector $\underline{b}$ are assured to converge if and only if all the eigenvalues of $\tilde{\boldsymbol{H}}$ have magnitude less than $1$, where $\tilde{\boldsymbol{H}}$ is defined in Lemma 6. When this condition is fulfilled, the mean values converge to the coordinates of the appropriate lattice point: $\underline{m}^{(\infty)}=\boldsymbol{G}\cdot\underline{b}$.

Proof:

Immediate from Lemma 6. ∎

Note that without adding random signs to the LDLC nonzero values, the all-ones vector will be an eigenvector of $\tilde{\boldsymbol{H}}$ with eigenvalue $\frac{\sum_{i=2}^{d}h_{i}}{h_{1}}$, which may exceed $1$.

Interestingly, recursion (6) is also obeyed by the error of the Jacobi method for solving systems of sparse linear equations [22] (see also Section VIII-A), when it is used to solve $\boldsymbol{H}\underline{m}=\underline{b}$ (with solution $\underline{m}=\boldsymbol{G}\underline{b}$). Therefore, the LDLC decoder can be viewed as a superposition of Jacobi solvers, one for each possible value of the integer valued vector $\underline{b}$.

We shall now turn to the convergence of the mean values of the wide messages. The asymptotic behavior is summarized in the following lemma.

Lemma 7

For a magic square LDLC with dimension $n$ and $\alpha<1$, consider only consistent Gaussians that relate to a given integer vector $\underline{b}$ and belong to wide messages. Denote the mean value of such a Gaussian that is sent from variable node $i$ at iteration $t$ by $m^{(t)}_{i}$, and arrange all these mean values in a column vector $\underline{m}^{(t)}$ of dimension $n$. Define the error vector $\underline{e}^{(t)}\stackrel{{\scriptstyle\Delta}}{{=}}\underline{m}^{(t)}-\boldsymbol{G}\underline{b}$. Then, for large $t$, $\underline{e}^{(t)}$ satisfies:

[TABLE]

where $\underline{y}$ is the noisy codeword and $\boldsymbol{F}$ is an $n\times n$ matrix defined by:

[TABLE]

Proof:

See Appendix E, where an alternative way to construct $\boldsymbol{F}$ from $\boldsymbol{H}$ is also presented. ∎

The conditions for convergence and steady state solution for the wide messages are described in the following theorem.

Theorem 3

For a magic square LDLC with $\alpha<1$, the mean values of the consistent Gaussians of the wide variable node messages that relate to a given integer vector $\underline{b}$ are assured to converge if and only if all the eigenvalues of $\boldsymbol{F}$ have magnitude less than $1$, where $\boldsymbol{F}$ is defined in Lemma 7. When this condition is fulfilled, the steady state solution is $\underline{m}^{(\infty)}=\boldsymbol{G}\cdot\underline{b}+(1-\alpha)(\boldsymbol{I}+\boldsymbol{F})^{-1}(\underline{y}-\boldsymbol{G}\cdot\underline{b})$.

Proof:

Immediate from Lemma 7. ∎

Unlike the narrow messages, the mean values of the wide messages do not converge to the appropriate lattice point coordinates. The steady state error depends on the difference between the noisy observation and the lattice point, as well as on $\alpha$, and it decreases to zero as $\alpha\rightarrow 1$. Note that the final PDF of a variable is generated by multiplying all the $d$ check node messages that arrive to the appropriate variable node. $d-1$ of these messages are wide, and therefore their mean values have a steady state error. One message is narrow, so it converges to an impulse at the lattice point coordinate. Therefore, the final product will be an impulse at the correct location, where the wide messages will only affect the magnitude of this impulse. As long as the mean values errors are not too large (relative to the width of the wide messages), this should not cause an impulse that corresponds to a wrong lattice point to have larger amplitude than the correct one. However, for large noise, these steady state errors may cause the decoder to deviate from the ML solution (As explained in Section IV-D).

To summarize the results for the mean values, we considered the mean values of all the consistent Gaussians that correspond to a given integer vector $\underline{b}$. A single Gaussian of this form exists in each of the $nd$ variable node messages that are generated in each iteration. For each variable node, $d-1$ messages are narrow (have variance that approaches zero) and a single message is wide (variance approaches a constant). Under certain conditions on $\boldsymbol{H}$, the mean values of all the narrow messages converge to the appropriate coordinate of the lattice point $\boldsymbol{G}\underline{b}$. Under additional conditions on $\boldsymbol{H}$, the mean values of the wide messages converge, but the steady state values contain an error term.

We analyzed the behavior of consistent Gaussian. It should be noted that there are many more non-consistent Gaussians. Furthermore non-consistent Gaussians are generated in each iteration for any existing consistent Gaussian. We conjecture that unless a Gaussian is consistent, or becomes consistent along the iterations, it fades out, at least at noise conditions where the algorithm converges. The reason is that non-consistency in the integer values leads to mismatch in the corresponding PDF’s, and so the amplitude of that Gaussian is attenuated.

We considered consistent Gaussians which correspond to a specific integer vector $\underline{b}$, but such a set of Gaussians exists for every possible choice of $\underline{b}$, i.e. for every lattice point. Therefore, the narrow messages will converge to a solution that has an impulse at the appropriate coordinate of every lattice point. This resembles the exact solution (1), so the key for proper convergence lies in the amplitudes: we would like the consistent Gaussians of the ML lattice point to have the largest amplitude for each message.

IV-D Convergence of the Amplitudes

We shall now analyze the behavior of the amplitudes of consistent Gaussians (as discussed later, this is not enough for complete convergence analysis, but it certainly gives insight to the nature of the convergence process and its properties). The behavior of the amplitudes of consistent Gaussians is summarized in the following lemma.

Lemma 8

For a magic square LDLC with dimension $n$, degree $d$ and $\alpha<1$, consider the $nd$ consistent Gaussians that relate to a given integer vector $\underline{b}$ in the variable node messages that are sent at iteration $t$ (one consistent Gaussian per message). Denote the amplitudes of these Gaussians by $p^{(t)}_{i}$, $i=1,2,...nd$, and define the log-amplitude as $l^{(t)}_{i}=\log p^{(t)}_{i}$. Arrange these $nd$ log-amplitudes in a column vector $\underline{l}^{(t)}$, such that element $(k-1)d+i$ corresponds to the message that is sent from variable node $k$ along an edge with weight $\pm h_{i}$. Assume further that the bipartite graph of the LDLC contains no 4-loops. Then, the log-amplitudes satisfy the following recursion:

[TABLE]

with initialization $\underline{l}^{(0)}=\underline{0}$. $\boldsymbol{A}$ is an $nd\times nd$ matrix which is all zeros except exactly $(d-1)^{2}$ ’$1$’s in each row and each column. The element of the excitation vector $\underline{a}^{(t)}$ at location $(k-1)d+i$ (where $k=1,2,...n$ and $i=1,2,...d$) equals:

[TABLE]

[TABLE]

where $\tilde{m}^{(t)}_{k,l}$ and $\tilde{V}^{(t)}_{k,l}$ denote the mean value and variance of the consistent Gaussian that relates to the integer vector $\underline{b}$ in the check node message that arrives to variable node $k$ at iteration $t$ along an edge with weight $\pm h_{l}$. $y_{k}$ is the noisy channel observation of variable node $k$, and $\hat{V}^{(t)}_{k,i}\stackrel{{\scriptstyle\Delta}}{{=}}\left(\frac{1}{\sigma^{2}}+\sum_{\begin{subarray}{c}l=1\\ l\neq i\end{subarray}}^{d}\frac{1}{\tilde{V}^{(t)}_{k,l}}\right)^{-1}$. Finally, $\underline{c}^{(t)}$ is a constant excitation term that is independent of the integer vector $\underline{b}$ (i.e. is the same for all consistent Gaussians). Note that an iteration is defined as sending variable node messages, followed by sending check node messages. The first iteration (where the variable nodes send the initialization PDF) is regarded as iteration [math].

Proof:

At the check node, the amplitude of a Gaussian at the convolution output is the product of the amplitudes of the corresponding Gaussians in the appropriate variable node messages. At the variable node, the amplitude of a Gaussian at the product output is the product of the amplitudes of the corresponding Gaussians in the appropriate check node messages, multiplied by the Gaussian scaling term, according to claim 2. Since we assume that the bipartite graph of the LDLC contains no 4-loops, an amplitude of a variable node message at iteration $t$ will therefore equal the product of $(d-1)^{2}$ amplitudes of Gaussians of variable node messages from iteration $t-1$, multiplied by the Gaussian scaling term. This proves (12) and shows that $A$ has $(d-1)^{2}$ ’1’s in every row. However, since each variable node message affects exactly $(d-1)^{2}$ variable node messages of the next iteration, $A$ must also have $(d-1)^{2}$ ’1’s in every column. The total excitation term $-\underline{a}^{(t)}-\underline{c}^{(t)}$ corresponds to the logarithm of the Gaussian scaling term. Each element of this scaling term results from the product of $d-1$ check node Gaussians and the channel Gaussian, according to claim 2. This scaling term sums over all the pairs of Gaussians, and in (13) the sum is separated to pairs that include the channel Gaussian and pairs that do not. The total excitation is divided between (13), which depends on the choice of the integer vector $\underline{b}$, and $\underline{c}^{(t)}$, which includes all the constant terms that are independent on $\underline{b}$ (including the normalization operation which is performed at the variable node). ∎

Since there are exactly $(d-1)^{2}$ ’$1$’s in each column of the matrix $\boldsymbol{A}$, it is easy to see that the all-ones vector is an eigenvector of $\boldsymbol{A}$, with eigenvalue $(d-1)^{2}$. If $d>2$, this eigenvalue is larger than $1$, meaning that the recursion (12) is non-stable.

It can be seen that the excitation term $\underline{a}^{(t)}$ has two components. The first term sums the squared differences between the mean values of all the possible pairs of received check node messages (weighted by the inverse product of the appropriate variances). It therefore measures the mismatch between the incoming messages. This mismatch will be small if the mean values of the consistent Gaussians converge to the coordinates of a lattice point (any lattice point). The second term sums the squared differences between the mean values of the incoming messages and the noisy channel output $y_{k}$. This term measures the mismatch between the incoming messages and the channel measurement. It will be smallest if the mean values of the consistent Gaussians converge to the coordinates of the ML lattice point.

The following lemma summarizes some properties of the excitation term $\underline{a}^{(t)}$.

Lemma 9

For a magic square LDLC with dimension $n$, degree $d$, $\alpha<1$ and no 4-loops, consider the consistent Gaussians that correspond to a given integer vector $\underline{b}$. According to Lemma 8, their amplitudes satisfy recursion (12). The excitation term $\underline{a}^{(t)}$ of (12), which is defined by (13), satisfies the following properties:

$a^{(t)}_{i}$, the $i$’th element of $\underline{a}^{(t)}$, is non-negative, finite and bounded for every $i$ and every $t$. Moreover, $a^{(t)}_{i}$ converges to a finite non-negative steady state value as $t\rightarrow\infty$. 2. 2.

$\lim_{t\rightarrow\infty}\sum_{i=1}^{nd}a^{(t)}_{i}=\frac{1}{2\sigma^{2}}(\boldsymbol{G}\underline{b}-\underline{y})^{T}\boldsymbol{W}(\boldsymbol{G}\underline{b}-\underline{y})$, where $\underline{y}$ is the noisy received codeword and $\boldsymbol{W}$ is a positive definite matrix defined by:

[TABLE]

[TABLE]

where $\boldsymbol{F}$ is defined in Lemma 7. 3. 3.

For an LDLC with degree $d>2$, the weighted infinite sum $\sum_{j=0}^{\infty}\frac{\sum_{i=1}^{nd}a^{(j)}_{i}}{(d-1)^{2j+2}}$ converges to a finite value.

Proof:

See Appendix F. ∎

The following theorem addresses the question of which consistent Gaussian will have the maximal asymptotic amplitude. We shall first consider the case of an LDLC with degree $d>2$, and then consider the special case of $d=2$ in a separate theorem.

Theorem 4

For a magic square LDLC with dimension $n$, degree $d>2$, $\alpha<1$ and no 4-loops, consider the $nd$ consistent Gaussians that relate to a given integer vector $\underline{b}$ in the variable node messages that are sent at iteration $t$ (one consistent Gaussian per message). Denote the amplitudes of these Gaussians by $p^{(t)}_{i}$, $i=1,2,...nd$, and define the product-of-amplitudes as $P^{(t)}\stackrel{{\scriptstyle\Delta}}{{=}}\prod_{i=1}^{nd}p^{(t)}_{i}$. Define further $S=\sum_{j=0}^{\infty}\frac{\sum_{i=1}^{nd}a^{(j)}_{i}}{(d-1)^{2j+2}}$, where $a^{(j)}_{i}$ is defined by (13) ($S$ is well defined according to Lemma 9). Then:

The integer vector $\underline{b}$ for which the consistent Gaussians will have the largest asymptotic product-of-amplitudes $\lim_{t\rightarrow\infty}P^{(t)}$ is the one for which $S$ is minimized. 2. 2.

The product-of-amplitudes for the consistent Gaussians that correspond to all other integer vectors will decay to zero in a super-exponential rate.

Proof:

As in Lemma 8, define the log-amplitudes $l^{(t)}_{i}\stackrel{{\scriptstyle\Delta}}{{=}}\log p^{(t)}_{i}$. Define further $s^{(t)}\stackrel{{\scriptstyle\Delta}}{{=}}\sum_{i=1}^{nd}l^{(t)}_{i}$. Taking the element-wise sum of (12), we get:

[TABLE]

with initialization $s^{(0)}=0$. Note that we ignored the term $\sum_{i=1}^{nd}c^{(t)}_{i}$. As shown below, we are looking for the vector $\underline{b}$ that maximizes $s^{(t)}$. Since (15) is a linear difference equation, and the term $\sum_{i=1}^{nd}c^{(t)}_{i}$ is independent of $\underline{b}$, its effect on $s^{(t)}$ is common to all $\underline{b}$ and is therefore not interesting.

Define now $\tilde{s}^{(t)}\stackrel{{\scriptstyle\Delta}}{{=}}\frac{s^{(t)}}{(d-1)^{2t}}$. Substituting in (15), we get:

[TABLE]

with initialization $\tilde{s}^{(0)}=0$, which can be solved to get:

[TABLE]

We would now like to compare the amplitudes of consistent Gaussians with various values of the corresponding integer vector $\underline{b}$ in order to find the lattice point whose consistent Gaussians will have largest product-of-amplitudes. From the definitions of $s^{(t)}$ and $\tilde{s}^{(t)}$ we then have:

[TABLE]

Consider two integer vectors $\underline{b}$ that relate to two lattice points. Denote the corresponding product-of-amplitudes by $P_{0}^{(t)}$ and $P_{1}^{(t)}$, respectively, and assume that for these two vectors $S$ converges to the values $S_{0}$ and $S_{1}$, respectively. Then, taking into account that $\lim_{t\rightarrow\infty}\tilde{s}^{(t)}=-S$, the asymptotic ratio of the product-of-amplitudes for these lattice points will be:

[TABLE]

It can be seen that if $S_{0}<S_{1}$, the ratio decreases to zero in a super exponential rate. This shows that the lattice point for which $S$ is minimized will have the largest product-of-amplitudes, where for all other lattice points, the product-of-amplitudes will decay to zero in a super-exponential rate (recall that the normalization operation at the variable node keeps the sum of all amplitudes in a message to be 1). This completes the proof of the theorem. ∎

We now have to find which integer valued vector $\underline{b}$ minimizes $S$. The analysis is difficult because the weighting factor inside the sum of (17) performs exponential weighting of the excitation terms, where the dominant terms are those of the first iterations. Therefore, we can not use the asymptotic results of Lemma 9, but have to analyze the transient behavior. However, the analysis is simpler for the case of an LDLC with row and column degree of $d=2$, so we shall first turn to this simple case (note that for this case, both the convolution in the check nodes and the product at the variable nodes involve only a single message).

Theorem 5

For a magic square LDLC with dimension $n$, degree $d=2$, $\alpha<1$ and no 4-loops, consider the $2n$ consistent Gaussians that relate to a given integer vector $\underline{b}$ in the variable node messages that are sent at iteration $t$ (one consistent Gaussian per message). Denote the amplitudes of these Gaussians by $p^{(t)}_{i}$, $i=1,2,...2n$, and define the product-of-amplitudes as $P^{(t)}\stackrel{{\scriptstyle\Delta}}{{=}}\prod_{i=1}^{2n}p^{(t)}_{i}$. Then:

The integer vector $\underline{b}$ for which the consistent Gaussians will have the largest asymptotic product-of-amplitudes $\lim_{t\rightarrow\infty}P^{(t)}$ is the one for which $(\boldsymbol{G}\underline{b}-\underline{y})^{T}\boldsymbol{W}(\boldsymbol{G}\underline{b}-\underline{y})$ is minimized, where $\boldsymbol{W}$ is defined by (14) and $\underline{y}$ is the noisy received codeword. 2. 2.

The product-of-amplitudes for the consistent Gaussians that correspond to all other integer vectors will decay to zero in an exponential rate.

Proof:

For $d=2$ (15) becomes:

[TABLE]

With solution:

[TABLE]

Denote $S_{a}=\lim_{j\rightarrow\infty}\sum_{i=1}^{2n}a^{(j)}_{i}$. $S_{a}$ is well defined according to Lemma 9. For large $t$, we then have $s^{(t)}\approx-t\cdot S_{a}$. Therefore, for two lattice points with excitation sum terms which approach $S_{a0},S_{a1}$, respectively, the ratio of the corresponding product-of-amplitudes will approach

[TABLE]

If $S_{a0}<S_{a1}$, the ratio decreases to zero exponentially (unlike the case of $d>2$ where the rate was super-exponential, as in (19)). This shows that the lattice point for which $S_{a}$ is minimized will have the largest product-of-amplitudes, where for all other lattice points, the product-of-amplitudes will decay to zero in an exponential rate (recall that the normalization operation at the variable node keeps the sum of all amplitudes in a message to be 1). This completes the proof of the second part of the theorem.

We still have to find the vector $\underline{b}$ that minimizes $S_{a}$. The basic difference between the case of $d=2$ and the case of $d>2$ is that for $d>2$ we need to analyze the transient behavior of the excitation terms, where for $d=2$ we only need to analyze the asymptotic behavior, which is much easier to handle.

According to Lemma 9, we have:

[TABLE]

where $\boldsymbol{W}$ is defined by (14) and $\underline{y}$ is the noisy received codeword. Therefore, for $d=2$, the lattice points whose consistent Gaussians will have largest product-of-amplitudes is the point for which $(\boldsymbol{G}\underline{b}-\underline{y})^{T}\boldsymbol{W}(\boldsymbol{G}\underline{b}-\underline{y})$ is minimized. This completes the proof of the theorem. ∎

For $d=2$ we could find an explicit expression for the “winning” lattice point. As discussed above, we could not find an explicit expression for $d>2$, since the result depends on the transient behavior of the excitation sum term, and not only on the steady state value. However, a reasonable conjecture is to assume that $\underline{b}$ that maximizes the steady state excitation will also maximize the term that depends on the transient behavior. This means that a reasonable conjecture is to assume that the “winning” lattice point for $d>2$ will also minimize an expression of the form (23).

Note that for $d>2$ we can still show that for “weak” noise, the ML point will have the minimal $S$. To see that, it comes out from (13) that for zero noise, the ML lattice point will have $a^{(t)}_{i}=0$ for every $t$ and $i$, where all other lattice points will have $a^{(t)}_{i}>0$ for at least some $i$ and $t$. Therefore, the ML point will have a minimal excitation term along the transient behavior so it will surely have the minimal $S$ and the best product-of-amplitudes. As the noise increases, it is difficult to analyze the transient behavior of $a^{(t)}_{i}$, as discussed above.

Note that the ML solution minimizes $(\boldsymbol{G}\underline{b}-\underline{y})^{T}(\boldsymbol{G}\underline{b}-\underline{y})$, where the above analysis yields minimization of $(\boldsymbol{G}\underline{b}-\underline{y})^{T}\boldsymbol{W}(\boldsymbol{G}\underline{b}-\underline{y})$. Obviously, for zero noise (i.e. $\underline{y}=\boldsymbol{G}\cdot\ \underline{b}$) both minimizations will give the correct solution with zero score. As the noise increases, the solutions may deviate from one another. Therefore, both minimizations will give the same solution for “weak” noise but may give different solutions for “strong” noise.

An example for another decoder that performs this form of minimization is the linear detector, which calculates $\hat{\underline{b}}=\left\lfloor\boldsymbol{H}\cdot\underline{y}\right\rceil$ (where $\left\lfloor x\right\rceil$ denotes the nearest integer to $x$). This is equivalent to minimizing $(\boldsymbol{G}\underline{b}-\underline{y})^{T}\boldsymbol{W}(\boldsymbol{G}\underline{b}-\underline{y})$ with $\boldsymbol{W}=\boldsymbol{H}^{T}\boldsymbol{H}=\boldsymbol{G}^{{-1}^{T}}\boldsymbol{G}^{-1}$. The linear detector fails to yield the ML solution if the noise is too strong, due to its inherent noise amplification.

For the LDLC iterative decoder, we would like that the deviation from the ML decoder due to the $\boldsymbol{W}$ matrix would be negligible in the expected range of noise variance. Experimental results (see Section IX) show that the iterative decoder indeed converges to the ML solution for noise variance values that approach channel capacity. However, for quantization or shaping applications (see Section VIII-B), where the effective noise is uniformly distributed along the Voronoi cell of the lattice (and is much stronger than the noise variance at channel capacity) the iterative decoder fails, and this can be explained by the influence of the $\boldsymbol{W}$ matrix on the minimization, as described above. Note from (14) that as $\alpha\rightarrow 1$, $\boldsymbol{W}$ approaches a scaled identity matrix, which means that the minimization criterion approaches the ML criterion. However, the variances converge as $\alpha^{t}$, so as $\alpha\rightarrow 1$ convergence time approaches infinity.

Until this point, we concentrated only on consistent Gaussians, and checked what lattice point maximizes the product-of-amplitudes of all the corresponding consistent Gaussians. However, this approach does not necessarily lead to the lattice point that will be finally chosen by the decoder, due to 3 main reasons:

It comes out experimentally that the strongest Gaussian in each message is not necessarily a consistent Gaussian, but a Gaussian that started as non-consistent and became consistent at a certain iteration. Such a Gaussian will finally converge to the appropriate lattice point, since the convergence of the mean values is independent of initial conditions. The non-consistency at the first several iterations, where the mean values are still very noisy, allows these Gaussians to accumulate stronger amplitudes than the consistent Gaussians (recall that the exponential weighting in (17) for $d>2$ results in strong dependency on the behavior at the first iterations). 2. 2.

There is an exponential number of Gaussians that start as non-consistent and become consistent (with the same integer vector $\underline{b}$) at a certain iteration, and the final amplitude of the Gaussians at the lattice point coordinates will be determined by the sum of all these Gaussians. 3. 3.

We ignored non-consistent Gaussians that endlessly remain non-consistent. We have not shown it analytically, but it is reasonable to assume that the excitation terms for such Gaussians will be weaker than for Gaussians that become consistent at some point, so their amplitude will fade away to zero. However, non-consistent Gaussians are born every iteration, even at steady state. The “newly-born” non-consistent Gaussians may appear as sidelobes to the main impulse, since it may take several iterations until they are attenuated. Proper choice of the coefficients of $\boldsymbol{H}$ may minimize this effect, as discussed in Sections III-A and V-A. However, these Gaussians may be a problem for small $d$ (e.g. $d=2$) where the product step at the variable node does not include enough messages to suppress them.

Note that the first two issues are not a problem for $d=2$, where the winning lattice point depends only on the asymptotic behavior. The amplitude of a sum of Gaussians that converged to the same coordinates will still be governed by (21) and the winning lattice point will still minimize (23). The third issue is a problem for small $d$, but less problematic for large $d$, as described above.

As a result, we can not regard the convergence analysis of the consistent Gaussians’ amplitudes as a complete convergence analysis. However, it can certainly be used as a qualitative analysis that gives certain insights to the convergence process. Two main observations are:

The narrow variable node messages tend to converge to single impulses at the coordinates of a single lattice point. This results from (19), (22), which show that the “non-winning” consistent Gaussians will have amplitudes that decrease to zero relative to the amplitude of the “winning” consistent Gaussian. This result remains valid for the sum of non-consistent Gaussians that became consistent at a certain point, because it results from the non-stable nature of the recursion (12), which makes strong Gaussians stronger in an exponential manner. The single impulse might be accompanied by weak “sidelobes” due to newly-born non-consistent Gaussians.

Interestingly, this form of solution is different from the exact solution (1), where every lattice point is represented by an impulse at the appropriate coordinate, with amplitude that depends on the Euclidean distance of the lattice point from the observation. The iterative decoder’s solution has a single impulse that corresponds to a single lattice point, where all other impulses have amplitudes that decay to zero. This should not be a problem, as long as the ML point is the remaining point (see discussion above). 2. 2.

We have shown that for $d=2$ the strongest consistent Gaussians relate to $\underline{b}$ that minimizes an expression of the form $(\boldsymbol{G}\underline{b}-\underline{y})^{T}\boldsymbol{W}(\boldsymbol{G}\underline{b}-\underline{y})$. We proposed a conjecture that this is also true for $d>2$. We can further widen the conjecture to say that the finally decoded $\underline{b}$ (and not only the $\underline{b}$ that relates to strongest consistent Gaussians) minimizes such an expression. Such a conjecture can explain why the iterative decoder works well for decoding near channel capacity, but fails for quantization or shaping, where the effective noise variance is much larger.

IV-E Summary of Convergence Results

To summarize the convergence analysis, it was first shown that the variable node messages are Gaussian mixtures. Therefore, it is sufficient to analyze the sequences of variances, mean values and relative amplitudes of the Gaussians in each mixture. Starting with the variances, it was shown that with proper choice of the magic square LDLC generating sequence, each variable node generates $d-1$ “narrow” messages, whose variance decreases exponentially to zero, and a single “wide” message, whose variance reaches a finite value. Consistent Gaussians were then defined as Gaussians that their generation process always involved the same integer at the same check node. Consistent Gaussians can then be related to an integer vector $\underline{b}$ or equivalently to the lattice point $\boldsymbol{G}\underline{b}$. It was then shown that under appropriate conditions on $\boldsymbol{H}$, the mean values of consistent Gaussians that belong to narrow messages converge to the coordinates of the appropriate lattice point. The mean values of wide messages also converge to these coordinates, but with a steady state error. Then, the amplitudes of consistent Gaussians were analyzed. For $d=2$ it was shown that the consistent Gaussians with maximal product-of-amplitudes (over all messages) are those that correspond to an integer vector $\underline{b}$ than minimizes $(\boldsymbol{G}\underline{b}-\underline{y})^{T}\boldsymbol{W}(\boldsymbol{G}\underline{b}-\underline{y})$, where $\boldsymbol{W}$ is a positive definite matrix that depends only on $\boldsymbol{H}$. The product-of-amplitudes for all other consistent Gaussians decays to zero. For $d>2$ the analysis is complex and depends on the transient behavior of the mean values and variances (and not only on their steady state values), but a reasonable conjecture is to assume that a same form of criterion is also minimized for $d>2$. The result is different from the ML lattice point, which minimizes $\left\|\boldsymbol{G}\cdot\underline{b}-\underline{y}\right\|^{2}$, where both criteria give the same point for weak noise but may give different solutions for strong noise. This may explain the experiments where the iterative decoder is successful in decoding the ML point for the AWGN channel near channel capacity, but fails in quantization or shaping applications where the effective noise is much stronger. These results also show that the iterative decoder converges to impulses at the coordinates of a single lattice point. It was then explained that analyzing the amplitudes of consistent Gaussians is not sufficient, so these results can not be regarded as a complete convergence analysis. However, the analysis gave a set of necessary conditions on $\boldsymbol{H}$, and also led to useful insights to the convergence process.

V Code Design

V-A Choosing the Generating Sequence

We shall concentrate on magic square LDLC, since they have inherent diversity of the nonzero elements in each row and column, which was shown above to be beneficial. It still remains to choose the LDLC generating sequence $h_{1},h_{2},...h_{d}$. Assume that the algorithm converged, and each PDF has a peak at the desired value. When the periodic functions are multiplied at a variable node, the correct peaks will then align. We would like that all the other peaks will be strongly attenuated, i.e. there will be no other point where the peaks align. This resembles the definition of the least common multiple (LCM) of integers: if the periods were integers, we would like to have their LCM as large as possible. This argument suggests the sequence $\{1/2,1/3,1/5,1/7,1/11,1/13,1/17,...\}$, i.e. the reciprocals of the smallest $d$ prime numbers. Since the periods are $1/h_{1},1/h_{2},...1/h_{d}$, we will get the desired property. Simulations have shown that increasing $d$ beyond $7$ with this choice gave negligible improvement. Also, performance was improved by adding some “dither” to the sequence, resulting in $\{1/2.31,1/3.17,1/5.11,1/7.33,1/11.71,1/13.11,1/17.55\}$. For $d<7$, the first $d$ elements are used.

An alternative approach is a sequence of the form $\{1,\epsilon,\epsilon,...,\epsilon\}$, where $\epsilon<<1$. For this case, every variable node will receive a single message with period $1$ and $d-1$ messages with period $1/\epsilon$. For small $\epsilon$, the period of these $d-1$ messages will be large and multiplication by the channel Gaussian will attenuate all the unwanted replicas. The single remaining replica will attenuate all the unwanted replicas of the message with period 1. A convenient choice is $\epsilon=\frac{1}{\sqrt{d}}$, which ensures that $\alpha=\frac{d-1}{d}<1$, as required by Theorem 1. As an example, for $d=7$ the sequence will be $\{1,\frac{1}{\sqrt{7}},\frac{1}{\sqrt{7}},\frac{1}{\sqrt{7}},\frac{1}{\sqrt{7}},\frac{1}{\sqrt{7}},\frac{1}{\sqrt{7}}\}$.

V-B Necessary Conditions on $\boldsymbol{H}$

The magic square LDLC definition and convergence analysis imply four necessary conditions on $\boldsymbol{H}$:

$|det(\boldsymbol{H})|=1$. This condition is part of the LDLC definition, which ensures proper density of the lattice points in $\mathbb{R}^{m}$. If $|det(\boldsymbol{H})|\neq 1$, it can be easily normalized by dividing $\boldsymbol{H}$ by $\sqrt[n]{|det(\boldsymbol{H})|}$. Note that practically we can allow $|det(\boldsymbol{H})|\neq 1$ as long as $\sqrt[n]{|det(\boldsymbol{H})|}\approx 1$, since $\sqrt[n]{|det(\boldsymbol{H})|}$ is the gain factor of the transmitted codeword. For example, if $n=1000$, having $|det(\boldsymbol{H})|=0.01$ is acceptable, since we have $\sqrt[n]{|det(\boldsymbol{H})|}=0.995$, which means that the codeword has to be further amplified by $20\cdot\log_{10}(0.995)=0.04$ dB, which is negligible.

Note that normalizing $\boldsymbol{H}$ is applicable only if $\boldsymbol{H}$ is non-singular. If $\boldsymbol{H}$ is singular, a row and a column should be sequentially omitted until $\boldsymbol{H}$ becomes full rank. This process may result in slightly reducing $n$ and a slightly different row and column degrees than originally planned. 2. 2.

$\alpha<1$, where $\alpha\stackrel{{\scriptstyle\Delta}}{{=}}\frac{\sum_{i=2}^{d}h_{i}^{2}}{h_{1}^{2}}$. This guarantees exponential convergence rate for the variances (Theorem 1). Choosing a smaller $\alpha$ results in faster convergence, but we should not take $\alpha$ too small since the steady state variance of the wide variable node messages, as well as the steady state error of the mean values of these messages, increases when $\alpha$ decreases, as discussed in Section IV-C. This may result in deviation of the decoded codeword from the ML codeword, as discussed in Section IV-D. For the first LDLC generating sequence of the previous subsection, we have $\alpha=0.92$ and $0.87$ for $d=7$ and $5$, respectively, which is a reasonable trade off. For the second sequence type we have $\alpha=\frac{d-1}{d}$. 3. 3.

All the eigenvalues of $\tilde{\boldsymbol{H}}$ must have magnitude less than $1$, where $\tilde{\boldsymbol{H}}$ is defined in Theorem 2. This is a necessary condition for convergence of the mean values of the narrow messages. Note that adding random signs to the nonzero $\boldsymbol{H}$ elements is essential to fulfill this necessary condition, as explained in Section IV-C. 4. 4.

All the eigenvalues of $\boldsymbol{F}$ must have magnitude less than $1$, where $\boldsymbol{F}$ is defined in Theorem 3. This is a necessary condition for convergence of the mean values of the wide messages.

Interestingly, it comes out experimentally that for large codeword length $n$ and relatively small degree $d$ (e.g. $n\geq 1000$ and $d\leq 10$), a magic square LDLC with generating sequence that satisfies $h_{1}=1$ and $\alpha<1$ results in $\boldsymbol{H}$ that satisfies all these four conditions: $\boldsymbol{H}$ is nonsingular without any need to omit rows and columns, $\sqrt[n]{|det(\boldsymbol{H})|}\approx 1$ without any need for normalization, and all eigenvalues of $\tilde{\boldsymbol{H}}$ and $\boldsymbol{F}$ have magnitude less than 1 (typically, the largest eigenvalue of $\tilde{\boldsymbol{H}}$ or $\boldsymbol{F}$ has magnitude of $0.94-0.97$, almost independently of $n$ and the choice of nonzero $\boldsymbol{H}$ locations). Therefore, by simply dividing the first generating sequence of the previous subsection by its first element, the constructed $\boldsymbol{H}$ meets all the necessary conditions, where the second type of sequence meets the conditions without any need for modifications.

V-C Construction of $\boldsymbol{H}$ for Magic Square LDLC

We shall now present a simple algorithm for constructing a parity check matrix for a magic square LDLC. If we look at the bipartite graph, each variable node and each check node has $d$ edges connected to it, one with every possible weight $h_{1},h_{2},...h_{d}$. All the edges that have the same weight $h_{j}$ form a permutation from the variable nodes to the check nodes (or vice versa). The proposed algorithm generates $d$ random permutations and then searches sequentially and cyclically for 2-loops (two parallel edges from a variable node to a check node) and 4-loops (two variable nodes that both are connected to a pair of check nodes). When such a loop is found, a pair is swapped in one of the permutations such that the loop is removed. A detailed pseudo-code for this algorithm is given in Appendix G.

VI Decoder Implementation

Each PDF should be approximated with a discrete vector with resolution $\Delta$ and finite range. According to the Gaussian Q-function, choosing a range of, say, $6\sigma$ to both sides of the noisy channel observation will ensure that the error probability due to PDF truncation will be $\approx 10^{-9}$. Near capacity, $\sigma^{2}\approx\frac{1}{2\pi e}$, so $12\sigma\approx 3$. Simulation showed that resolution errors became negligible for $\Delta=1/64$. Each PDF was then stored in a $L=256$ elements vector, corresponding to a range of size $4$.

At the check node, the PDF $f_{j}(x)$ that arrives from variable node $j$ is first expanded by $h_{j}$ (the appropriate coefficient of $\boldsymbol{H}$) to get $f_{j}(x/h_{j})$. In a discrete implementation with resolution $\Delta$ the PDF is a vector of values $f_{j}(k\Delta)$, $k\in\mathbb{Z}$. As described in Section V, we shall usually use $h_{j}\leq 1$ so the expanded PDF will be shorter than the original PDF. If the expand factor $1/|h_{j}|$ was an integer, we could simply decimate $f_{j}(k\Delta)$ by $1/|h_{j}|$. However, in general it is not an integer so we should use some kind of interpolation. The PDF $f_{j}(x)$ is certainly not band limited, and as the iterations go on it approaches an impulse, so simple interpolation methods (e.g. linear) are not suitable. Suppose that we need to calculate $f_{j}((k+\epsilon)\Delta)$, where $-0.5\leq\epsilon\leq 0.5$. A simple interpolation method which showed to be effective is to average $f_{j}(x)$ around the desired point, where the averaging window length $l_{w}$ is chosen to ensure that every sample of $f_{j}(x)$ is used in the interpolation of at least one output point. This ensures that an impulse can not be missed. The interpolation result is then $\frac{1}{2l_{w}+1}\sum_{i=-l_{w}}^{l_{w}}f_{j}((k-i)\Delta)$, where $l_{w}=\left\lfloor\frac{\left\lceil 1/|h_{j}|\right\rceil}{2}\right\rfloor$.

The most computationally extensive step at the check nodes is the calculation the convolution of $d-1$ expanded PDF’s. An efficient method is to calculate the fast Fourier transforms (FFTs) of all the PDF’s, multiply the results and then perform inverse FFT (IFFT). The resolution of the FFT should be larger than the expected convolution length, which is roughly $L_{out}\approx L\cdot\sum_{i=1}^{d}h_{i}$, where $L$ denotes the original PDF length. Appendix H shows a way to use FFTs of size $1/\Delta$, where $\Delta$ is the resolution of the PDF. Usually $1/\Delta<<L_{out}$ so FFT complexity is significantly reduced. Practical values are $L=256$ and $\Delta=1/64$, which give an improvement factor of at least $4$ in complexity.

Each variable node receives $d$ check node messages. The output variable node message is calculated by generating the product of $d-1$ input messages and the channel Gaussian. As the iterations go on, the messages get narrow and may become impulses, with only a single nonzero sample. Quantization effects may cause impulses in two messages to be shifted by one sample. This will result in a zero output (instead of an impulse). Therefore, it was found useful to widen each check node message $Q(k)$ prior to multiplication, such that $Q_{w}(k)=\sum_{i=-1}^{1}Q(k+i)$, i.e. the message is added to its right shifted and left shifted (by one sample) versions.

VII Computational Complexity and Storage Requirements

Most of the computational effort is invested in the $d$ FFT’s and $d$ IFFT’s (of length $1/\Delta$) that each check node performs each iteration. The total number of multiplications for $t$ iterations is $o\left(n\cdot d\cdot t\cdot\frac{1}{\Delta}\cdot\log_{2}(\frac{1}{\Delta})\right)$. As in binary LDPC codes, the computational complexity has the attractive property of being linear with block length. However, the constant that precedes the linear term is significantly higher, mainly due to the FFT operations.

The memory requirements are governed by the storage of the $nd$ check node and variable node messages, with total memory of $o(n\cdot d\cdot L)$. Compared to binary LDPC, the factor of $L$ significantly increases the required memory. For example, for $n=10,000$, $d=7$ and $L=256$, the number of storage elements is of the order of $10^{7}$.

VIII Encoding and Shaping

VIII-A Encoding

The LDLC encoder has to calculate $\underline{x}=\boldsymbol{G}\cdot\underline{b}$, where $\underline{b}$ is an integer message vector. Note that unlike $\boldsymbol{H}$, $\boldsymbol{G}=\boldsymbol{H}^{-1}$ is not sparse, in general, so the calculation requires computational complexity and storage of $o(n^{2})$. This is not a desirable property because the decoder’s computational complexity is only $o(n)$. A possible solution is to use the Jacobi method [22] to solve $\boldsymbol{H}\cdot\underline{x}=\underline{b}$, which is a system of sparse linear equations. Using this method, a magic square LDLC encoder calculates several iterations of the form:

[TABLE]

with initialization $\underline{x}^{(0)}=\underline{0}$. The matrix $\tilde{\boldsymbol{H}}$ is defined in Lemma 6 of Section IV-C. The vector $\underline{\tilde{b}}$ is a permuted and scaled version of the integer vector $\underline{b}$, such that the $i$’th element of $\underline{\tilde{b}}$ equals the element of $\underline{b}$ for which the appropriate row of $\boldsymbol{H}$ has its largest magnitude value at the $i$’th location. This element is further divided by this largest magnitude element.

A necessary and sufficient condition for convergence to $\underline{x}=\boldsymbol{G}\cdot\underline{b}$ is that all the eigenvalues of $\tilde{\boldsymbol{H}}$ have magnitude less than $1$ [22]. However, it was shown that this is also a necessary condition for convergence of the LDLC iterative decoder (see Sections IV-C, V-B), so it is guaranteed to be fulfilled for a properly designed magic square LDLC. Since $\tilde{\boldsymbol{H}}$ is sparse, this is an $o(n)$ algorithm, both in complexity and storage.

VIII-B Shaping

For practical use with the power constrained AWGN channel, the encoding operation must be accompanied by shaping, in order to prevent the transmitted codeword’s power from being too large. Therefore, instead of mapping the information vector $\underline{b}$ to the lattice point $\underline{x}=\boldsymbol{G}\cdot\underline{b}$, it should be mapped to some other lattice point $\underline{x}^{\prime}=\boldsymbol{G}\cdot\underline{b}^{\prime}$, such that the lattice points that are used as codewords belong to a shaping region (e.g. an $n$-dimensional sphere). The shaping operation is the mapping of the integer vector $\underline{b}$ to the integer vector $\underline{b}^{\prime}$.

As explained in Section II-A, this work concentrates on the lattice design and the lattice decoding algorithm, and not on the shaping region or shaping algorithms. Therefore, this section will only highlight some basic shaping principles and ideas.

A natural shaping scheme for lattice codes is nested lattice coding [12]. In this scheme, shaping is done by quantizing the lattice point $\boldsymbol{G}\cdot\underline{b}$ onto a coarse lattice $\boldsymbol{G}^{\prime}$, where the transmitted codeword is the quantization error, which is uniformly distributed along the Voronoi cell of the coarse lattice. If the second moment of this Voronoi cell is close to that of an $n$-dimensional sphere, the scheme will attain close-to-optimal shaping gain. Specifically, assume that the information vector $\underline{b}$ assumes integer values in the range $0,1,...M-1$ for some constant integer $M$. Then, we can choose the coarse lattice to be $\boldsymbol{G^{\prime}}=M\boldsymbol{G}$. The volume of the Voronoi cell for this lattice is $M^{n}$, since we assume $det(\boldsymbol{G})=1$ (see Section II-A). If the shape of the Voronoi cell resembles an $n$-dimensional sphere (as expected from a capacity approaching lattice code), it will attain optimal shaping gain (compared to uncoded transmission of the original integer sequence $\underline{b}$).

The shaping operation will find the coarse lattice point $M\boldsymbol{G}\underline{k}$, $\underline{k}\in\mathbb{Z}^{n}$, which is closest to the fine lattice point $\underline{x}=\boldsymbol{G}\cdot\underline{b}$. The transmitted codeword will be:

[TABLE]

where $\underline{b}^{\prime}\stackrel{{\scriptstyle\Delta}}{{=}}\underline{b}-M\underline{k}$ (note that the “inverse shaping” at the decoder, i.e. transforming from $\underline{b}^{\prime}$ to $\underline{b}$, is a simple modulo calculation: $\underline{b}=\underline{b}^{\prime}\mod M$). Finding the closest coarse lattice point $M\boldsymbol{G}\underline{k}$ to $\underline{x}$ is equivalent to finding the closest fine lattice point $\boldsymbol{G}\cdot\underline{k}$ to the vector $\underline{x}/M$. This is exactly the operation of the iterative LDLC decoder, so we could expect that is could be used for shaping. However, simulations show that the iterative decoding finds a vector $\underline{k}$ with poor shaping gain. The reason is that for shaping, the effective noise is much stronger than for decoding, and the iterative decoder fails to find the nearest lattice point if the noise is too large (see Section IV-D).

Therefore, an alternative algorithm has to be used for finding the nearest coarse lattice point. The complexity of finding the nearest lattice point grows exponentially with the lattice dimension $n$ and is not feasible for large dimensions [23]. However, unlike decoding, for shaping applications it is not critical to find the exact nearest lattice point, and approximate algorithms may be considered (see [15]). A possible method [24] is to perform QR decomposition on $\boldsymbol{G}$ in order to transform to a lattice with upper triangular generator matrix, and then use sequential decoding algorithms (such as the Fano algorithm) to search the resulting tree. The main disadvantage of this approach is computational complexity and storage of at least $o(n^{2})$. Finding an efficient shaping scheme for LDLC is certainly a topic for further research.

IX Simulation Results

Magic square LDLC with the first generating sequence of Section V-A (i.e. $\{1/2.31,1/3.17,1/5.11,1/7.33,1/11.71,1/13.11,1/17.55\}$) were simulated for the AWGN channel at various block lengths. The degree was $d=5$ for $n=100$ and $d=7$ for all other $n$. For $n=100$ the matrix $\boldsymbol{H}$ was further normalized to get $\sqrt[n]{det(\boldsymbol{H})}=1$. For all other $n$, normalizing the generating sequence such that the largest element has magnitude $1$ also gave the desired determinant normalization (see Section V-B). The $\boldsymbol{H}$ matrices were generated using the algorithm of Section V-C. PDF resolution was set to $\Delta=1/256$ with a total range of $4$, i.e. each PDF was represented by a vector of $L=1024$ elements. High resolution was used since our main target is to prove the LDLC concept and eliminate degradation due to implementation considerations. For this reason, the decoder was used with 200 iterations (though most of the time, a much smaller number was sufficient).

In all simulations the all-zero codeword was used. Approaching channel capacity is equivalent to $\sigma^{2}\rightarrow\frac{1}{2\pi e}$ (see Section II-A), so performance is measured in symbol error rate (SER), vs. the distance of the noise variance $\sigma^{2}$ from capacity (in dB). The results are shown in Figure 4. At SER of $10^{-5}$, for $n=100000$, $10000$, $1000$, $100$ we can work as close as $0.6$dB, $0.8$dB, $1.5$dB and $3.7$dB from capacity, respectively.

Similar results were obtained for $d=7$ with the second type of generating sequence of Section V-A, i.e. $\{1,\frac{1}{\sqrt{7}},\frac{1}{\sqrt{7}},\frac{1}{\sqrt{7}},\frac{1}{\sqrt{7}},\frac{1}{\sqrt{7}},\frac{1}{\sqrt{7}}\}$. Results were slightly worse than for the first generating sequence (by less than 0.1 dB). Increasing $d$ did not give any visible improvement.

X Conclusion

Low density lattice codes (LDLC) were introduced. LDLC are novel lattice codes that can approach capacity and be decoded efficiently. Good error performance within $\sim 0.5$dB from capacity at block length of 100,000 symbols was demonstrated. Convergence analysis was presented for the iterative decoder, which is not complete, but yields necessary conditions on $\boldsymbol{H}$ and significant insight to the convergence process. Code parameters were chosen from intuitive arguments, so it is reasonable to assume that when the code structure will be more understood, better parameters could be found, and channel capacity could be approached even closer.

Multi-input, multi-output (MIMO) communication systems have become popular in recent years. Lattice codes have been proposed in this context as space-time codes (LAST) [25]. The concatenation of the lattice encoder and the MIMO channel generates a lattice. If LDLC are used as lattice codes and the MIMO configuration is small, the inverse generator matrix of this concatenated lattice can be assumed to be sparse. Therefore, the MIMO channel and the LDLC can be jointly decoded using an LDLC-like decoder. However, even if a magic square LDLC is used as the lattice code, the concatenated lattice is not guaranteed to be equivalent to a magic square LDLC, and the necessary conditions for convergence are not guaranteed to be fulfilled. Therefore, the usage of LDLC for MIMO systems is a topic for further research.

Appendix A Exact PDF Calculations

Given the $n$-dimensional noisy observation $\underline{y}=\underline{x}+\underline{w}$ of the transmitted codeword $\underline{x}=\boldsymbol{G}\underline{b}$, we would like to calculate the probability density function (PDF) $f_{x_{k}|\underline{y}}(x_{k}|\underline{y})$. We shall start by calculating $f_{\underline{x}|\underline{y}}(\underline{x}|\underline{y})=\frac{f_{\underline{x}}(\underline{x})f_{\underline{y}|\underline{x}}(\underline{y}|\underline{x})}{f_{\underline{y}}(\underline{y})}$. Denote the shaping region by $B$ ($\boldsymbol{G}$ will be used to denote both the lattice and its generator matrix). $f_{\underline{x}}(\underline{x})$ is a sum of $|\boldsymbol{G}\cap B|$ $n$-dimensional Dirac delta functions, since $\underline{x}$ has nonzero probability only for the lattice points that lie inside the shaping region. Assuming further that all codewords are used with equal probability, all these delta functions have equal weight of $\frac{1}{|\boldsymbol{G}\cap B|}$. The expression for $f_{\underline{y}|\underline{x}}(\underline{y}|\underline{x})$