Gradient Flow Decoding for LDPC Codes

The code length is denoted by $n$. A binary sparse matrix over $\mathbb{F}_{2}$, $\bm{H}\in\mathbb{F}_{2}^{m\times n}$, is a parity check matrix for an LDPC code. The binary linear code $\tilde{C}(\bm{H})$ is defined by

A binary to bipolar transform $\beta:\mathbb{F}_{2}\rightarrow\{1,-1\}$ defined as $\beta(0)\equiv 1$ and $\beta(1)\equiv-1$ transforms $\tilde{C}(\bm{H})$ into the bipolar code defined by

The index sets $A(i)$ and $B(j)$ are defined as

respectively, where $H_{i,j}$ denotes the $(i,j)$-element of $\bm{H}$. The notation $[n]$ denotes the set of consecutive integers $\{1,2,\ldots,n\}$. The multivariate Gaussian distribution with mean vector $\bm{m}$ and covariance $\bm{\Sigma}$ is denoted by ${\cal N}(\bm{m},\bm{\Sigma})$.

The focus of current research lies in optimization-based approaches. Several works have been developed in this category, with the celebrated work by Feldman [9] on linear programming decoding providing a clear demonstration that the decoding problem can be viewed as an optimization problem.

Applications of interior point methods for solving a convex problem to decoding problems have been studied [10] [11]. A gradient descent formulation of a non-convex objective function leads to the GDBF algorithm [2].

Some variants of the GDBF algorithm, such as the noisy GDBF algorithm [12], have shown considerable improvement in decoding performance. Additionally, recent works have presented ADMM-based decoding algorithms for LDPC codes, including those by Zhang, Liu, and Wang [13][14][15][16], which have shown promising results. Although the advancement of these optimization-based decoding algorithms is significant, all of them are discrete-time algorithms.

Research on continuous-time dynamical systems for signal processing is a relatively new field and still in its early stages. However, it has shown promise in solving optimization problems and signal recovery tasks. For example, a MIMO MMSE detector based on continuous-time dynamics was proposed in [17], and a sparse signal recovery dynamics was studied in [18]. Further studies in this field, including the proposed work on LDPC decoding using gradient flow dynamics, can help to clarify the benefits and limitations of signal processing based on continuous dynamics.

To define the gradient flow, we need an appropriate potential energy function. The code potential energy function for $C(\bm{H})$ is a multivariate polynomial defined as

where $\bm{x}=(x_{1},\ldots,x_{n})^{T}\in\mathbb{R}^{n}$. The parameters $\alpha\in\mathbb{R}_{+}$ and $\beta\in\mathbb{R}_{+}$ control the strength of the constraints. In the right-hand side of this equation, the first term represents the bipolar constraint for $\bm{x}\in\{+1,-1\}^{n}$, and the second term corresponds to the parity constraint induced by $\bm{H}$, i.e., if $\bm{x}\in C(\bm{H})$, we have $\left(\prod_{j\in A(i)}x_{j}\right)-1=0$ for any $i\in[m]$. The code potential energy function was first introduced in the work on proximal decoding algorithm [19].

Since the polynomial $h(\bm{x})$ has a sum-of-squares form, the polynomial can be regarded as a penalty function that yields a positive penalty value for non-codeword vectors in $\mathbb{R}^{n}$. The code potential energy $h_{\alpha,\beta}(\bm{x})$ is inspired by the non-convex parity constraint function used in the GDBF objective function [2]. The sum-of-squares form directly implies the most important property of $h_{\alpha,\beta}(\bm{x})$, i.e., the inequality $h_{\alpha,\beta}(\bm{x})\geq 0$ holds for any $\bm{x}\in\mathbb{R}^{n}$. The equality holds if and only if $\bm{x}\in C(\bm{H})$.

In the following discussion, we need the gradient of $h_{\alpha,\beta}(\bm{x})$. The first-order derivative of $h_{\alpha,\beta}(\bm{x})$ with respect to $x_{k}(k\in[n])$ is given by

The gradient $\nabla h_{\alpha,\beta}(\bm{x})$ is thus given by

The point $\bm{z}\in\mathbb{R}^{n}$ satisfying the equality $\nabla h_{\alpha,\beta}(\bm{z})=\bm{0}$ is a stationary point of $h_{\alpha,\beta}$. For any codeword $\bm{x}\in C(\bm{H})$, $x_{k}^{2}=1$ for any $k\in[n]$ and $\prod_{j\in A(i)}x_{j}=1$ holds for any $i\in[m]$. This means $\nabla h_{\alpha,\beta}(\bm{x})=\bm{0}$, i.e., a codeword vector is a stationary point of $h$.

The gradient flow dynamics, also known as steepest descent dynamics, is a continuous-time dynamics defined by an ODE:

where $F:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is a potential energy function. The system’s state $\bm{x}:\mathbb{R}\rightarrow\mathbb{R}^{n}$ evolves to reduce the potential energy $F$ as the time variable $t$ increases. This continuous dynamical system can be viewed as the continuous counterpart of the gradient descent method. If the potential function $F$ is strictly convex, the equilibrium point of the dynamics coincides with the minimum point of $F$. Therefore, the gradient flow dynamics can be seen as a continuous-time minimization process of the potential energy. If $F$ is a non-convex function, then the gradient flow dynamics finds a stationary point of $F$ as $t\rightarrow\infty$.

We assume an AWGN channel in this paper. A transmitter send a codeword $\bm{s}\in C(\bm{H})$ and a receiver obtain a received word

where $\bm{n}\sim{\cal N}(\bm{0},\sigma^{2}\bm{I})$.

Let $f$ be a potential energy function defined by

The first term of the potential energy function can be regarded as the negative log likelihood function for an AWGN channel. The second term is a penalty term attracting $\bm{x}$ to a codeword of $C(\bm{H})$.

The gradient flow decoding is simply a gradient flow dynamics based on the potential energy function $f$ that is given by ${d\bm{x}(t)}/{dt}=-\nabla f(\bm{x}(t)).$ This ODE can be rewritten as follows.

In the above ODE, the notation $\bm{x}(t)$ is abbreviated to $\bm{x}$ for simplicity. If we have a suitably powerful analog circuit or analog computer, we can simulate the dynamics of this ODE. This simulation process is equivalent to a decoding process, where the equilibrium point represents the decoding result. A possible choice for the initial point is to set $\bm{x}_{0}=\bm{0}$. Alternatively, for certain cases, a small-norm vector or a random vector could be used as the initial point.

To evaluate the decoding performance of the GF decoding (11), a numerical method is required to solve the ODE. The simplest numerical method for solving simultaneous nonlinear differential equations is the Euler method [20]. While the convergence order of the Euler method is lower than that of higher-order methods, it is straightforward to use and can provide sufficiently accurate solutions if the time interval is sufficiently fine discretized. Therefore, in the present study, we will use the Euler method to solve the ODE (11).

Suppose that we require to simulate the dynamics defined by the above ODE (11) in the time interval $0\leq t\leq T$. The interval is divided into $N$ bins where $N$ denotes the number of discretized intervals. The discrete-time ticks $t_{k}=k\eta(k=0,1,\ldots,N)$ defines the boundaries of the bins where the width of a bin is given by $\eta\equiv T/N.$ It should be noted that the choice of the bin width $\eta$ is crucial in order to ensure the stability and the accuracy of the Euler method. By using the Euler method, the ODE (11) can be approximated by

for $k=0,1,2,\ldots,N-1$. The initial condition of this recursion is $\bm{x}^{(0)}=\bm{x}_{0}$.

Here, we present a small example to illustrate a decoding process. Suppose that we have a repetition code of length 2, i.e., $C={(+1,+1),(-1,-1)}$. The code potential energy function for $C$ is given by:

where the parameters $\alpha$ and $\beta$ are set to 1.

Assume that a transmitted word is $\bm{x}=(1,1)$ and that the corresponding received word is $\bm{y}=(0.6027,0.8244)$. In this case, the ODE (11) for the GF decoding becomes

The initial condition for the ODE is set to $\bm{x}(0)=(0,0)$. Figure 1 shows the solution curve of the ODE, with the white small circle representing the initial point and the black small circle indicating the equilibrium point $(0.9642,0.9901)$. The solution curve is obtained numerically using the Euler method. In Fig. 1, the contour curves of the potential energy function $f(\bm{x})$ are also plotted. We can observe that the solution curve is orthogonal to the contour curves because the solution path follows the negative gradient vector field of the potential energy. This means that the curve is a steepest descent curve for $f$. By rounding the equilibrium point, we obtain $\hat{\bm{x}}=(1,1)$, which is the correct estimated word.

If the state follows a gradient flow dynamics, the potential energy should decrease as the state evolves. To confirm the evolution of potential energy, we conducted an experiment as follows. For each trial, a random bipolar codeword $\bm{s}$ of a (3,6)-regular LDPC code was generated. The received word $\bm{y}$ was also randomly sampled as $\bm{y}=\bm{s}+\bm{n}$, where $\bm{n}\sim{\cal N}(\bm{0},\sigma^{2}\bm{I})$. The received vector $\bm{y}$ was then provided to the gradient flow decoder. We conducted 50 trials in total.

Figure 2 plots the potential energy $f(\bm{x}(t))$ as a function of time $t$. As expected, we observed that the value of potential energy rapidly decreased. In the regime $t\geq 1$, the decrement of the potential energy was invisible, likely because the state vector was sufficiently close to the equilibrium point in such a regime.

It may be informative to observe a decoding process of the GF decoding. Figure 3 presents a decoding process of the (3,6)-regular LDPC code. The dashed line represents the transmitted word where the horizontal axis represents the index of a vector from $1$ to $n=204$. The solid line in Fig. 3 indicates the value of $x_{i}$ in $\bm{x}$. The upper panel corresponds to time $t=0.1$ and lower panel corresponds to time $t=1$. At the beginning of the decoding process (upper panel), the magnitude of $x_{i}$ in the state vector $\bm{x}$ are fairly small. They are not sufficient to yield the final estimation. When $t=1$ (lower panel), the values of $x_{i}$ are significantly close to the original vector. The state vector may be close to the equilibrium point. This decoding result corresponds to successful decoding.

Figure 4 displays the solution curves of each element in the state vector $\bm{x}$, where $x_{i}(t)$ (for $i\in[204]$) represents the solution of the ODE (11) plotted as a function of time $t$. Notably, all solution curves are smooth and continuous. The left panel illustrates the solution curves for ${\sf SNR}=4$dB, where the equilibrium point is in proximity to a bipolar vector. Conversely, in the right panel (${\sf SNR}=0$dB), the solution curves are more widely dispersed.