Thus, it is sufficient for us to bound
$\sum_{i=1}^{T}\|x_{t}\|_{*}^{2}=\sum_{i=1}^{T}\phi_{i}$. Similar to Lemma 17
we have the following result

######  Lemma 19.

[27] When $M\leq\lambda$, we have

$\displaystyle\mathbb{E}[\phi_{i}]\leq(\sigma+\frac{M\sigma}{\lambda})^{2}\leq
4\sigma^{2},$ $\displaystyle\phi_{i}\leq(2M+\lambda)^{2}\leq 9\lambda^{2},$
$\displaystyle[\mathbb{E}(\phi_{i}^{2})]^{\frac{1}{2}}\leq(\sigma+\frac{M\sigma}{\lambda})(2M+\lambda)\leq
6\lambda\sigma.$

Thus, by Berstern’s inequality, we have if
$\tau=O\left(\frac{\sigma^{2}T}{\lambda^{2}}\right)$

$\displaystyle\text{Pr}[\sum_{t=1}^{T}||x_{t}||_{*}^{2}\geq
4\sigma^{2}T+18\lambda\sigma\sqrt{T\tau}]\leq\exp(-\frac{9\tau}{2+\frac{3\sqrt{\tau}\lambda}{\sigma\sqrt{T}}})\leq\exp(-\tau).$

In total, let $\gamma_{t}=\bar{\gamma}$, we have with probability at least
$1-2\exp(-\tau)$

$\mathbb{E}[\mathcal{L}(\hat{w})]-L(\theta^{*})\leq
O\left(\frac{D_{\Phi}(\theta^{*},w_{0})\cdot\bar{\gamma}}{T}+\frac{M\sigma^{2}}{\lambda}+\frac{\sigma
M\sqrt{\tau}}{\sqrt{T}}+\frac{\sigma^{2}}{\bar{\gamma}}+\frac{M\sigma\sqrt{\tau}}{\sqrt{T}\bar{\gamma}}+\frac{\log(\frac{n}{\delta})\cdot
dT\kappa(\beta
M+\lambda)^{2}\cdot\log(1/\delta)}{n^{2}\epsilon^{2}\bar{\gamma}}\right).$
(27)

Let $\frac{\bar{\gamma}}{T}=O(\frac{(\beta
M+\lambda)\sqrt{d\log(1/\delta)}}{nM\epsilon})$, and since
$D_{\Phi}(\theta^{*},w_{0})=\Phi(\theta^{*})\leq\frac{\kappa M^{2}}{2}$ we
have

$\mathbb{E}[\mathcal{L}(\hat{w})]-L(\theta^{*})\leq\tilde{O}\left(\frac{M\sigma^{2}}{\lambda}+\frac{\sigma
M\sqrt{\tau}}{\sqrt{T}}+\frac{M\sigma^{2}}{\bar{\gamma}}+\frac{(\beta
M+\lambda)M\kappa\sqrt{d\log(1/\delta)}}{n\epsilon}\right).$

Let
$\lambda=\frac{\sigma\sqrt{n\epsilon}}{\sqrt[4]{\kappa^{2}d\log(1/\delta)}}\geq\max\\{\beta,1\\}M$,
we have

$\mathbb{E}[\mathcal{L}(\hat{w})]-\mathcal{L}(\theta^{*})\leq
O\left(\frac{M\sigma\kappa\sqrt[4]{d\log(1/\delta)}}{\sqrt{n\epsilon}}+\frac{\sigma
M\sqrt{\tau}}{\sqrt{T}}+\frac{M\sigma^{2}}{\bar{\gamma}}\right).$

Let $\bar{\gamma}=\sqrt{T}$, then $\sqrt{T}=O(\frac{Mn\epsilon}{(\beta
M+\lambda)\sqrt{d\log(1/\delta)}})$, and it holds that

$\mathbb{E}[\mathcal{L}(\hat{w})]-\mathcal{L}(\theta^{*})\leq
O\left(\frac{M\max\\{\sigma^{2},\sigma\\}\sqrt[4]{\kappa^{2}d\log(1/\delta)}\sqrt{\log(1/\delta^{{}^{\prime}}})}{\sqrt{n\epsilon}}\right)$

w.p. at least $1-\delta^{{}^{\prime}}$. ∎

###  B.6 Proof of Theorem 11

Algorithm 7 Truncated DP Batched Mirror Descent

1:Input: Dataset $D$, loss function $\ell$, initial point $w_{0}=0$, smooth
parameter $\beta$ and and $\lambda$.

2:Divide the permuted data into $T$ batches $\\{B_{i}\\}_{i=1}^{T}$ where
$|B_{i}|=\frac{n}{T}$ for all $i=1,\cdots,T$

3:for $t=1,\cdots,T$ do

4: for each $x\in B_{t}$ do

5: $g_{x}=\begin{cases}\nabla\ell(w_{t-1},x)&\text{if}\leavevmode\nobreak\
||\nabla\ell(w_{t-1},x)||_{*}\leq\beta M+\lambda\\\
0&{\text{otherwise}}\end{cases}$

6: end for

7: Let

8:
$w_{t}=\underset{w\in\mathcal{\mathcal{C}}}{\arg\min}\left\\{\langle\frac{\underset{{x\in
B_{t}}}{\sum}g_{x}}{|B_{t}|}+Z^{t},w\rangle+\gamma_{t}\cdot
D_{\Phi}(w,w_{t-1})\right\\},$ where
$Z^{t}\sim\mathcal{GG}_{||\cdot||_{+}}(\sigma_{1}^{2})$ with
$\sigma_{1}^{2}=O\left(\frac{\kappa(\beta
M+\lambda)^{2}\cdot\log(1/\delta)}{|B_{t}|^{2}\epsilon^{2}}\right)$,
$||\cdot||_{+}$ is the smooth norm for $(\mathbf{E},||\cdot||_{*})$.
$\kappa=\min\\{\frac{1}{p-1},\log d\\}$ and
$\Phi(x)=\frac{\kappa}{2}||x||_{\kappa_{+}}^{2}$ with
$\kappa_{+}=\frac{\kappa}{\kappa-1}$.

9:end for

10:return
$\hat{w}=(\sum_{t=1}^{T}\gamma_{t}^{-1})^{-1}\cdot\sum_{t=1}^{T}\gamma_{t}^{-1}w_{t}$

We propose our method in Algorithm 7. Note that there are two key differences
compared to Algorithm 4. First, since we do not need the privacy amplification
via shuffling, there is no shuffling step. Secondly, instead of adding noise
to each truncated gradient $g_{x}$, here we add a generalized Gaussian noise
to the averages of the gradients for each batch. In the following we will
provide our theoretical results.

######  Theorem 16.

For the $\ell_{p}^{d}$ space with $1<p<2$, suppose Assumption 4 holds and
assume $n$ is large enough such that
$O((\frac{\sqrt{n\epsilon}M}{\kappa\sqrt[4]{d\log(1/\delta)}})^{\frac{2}{3}})\geq\max\\{\beta,1\\}M$.
For any $0<\epsilon,\delta<1$, Algorithm 7 is $(\epsilon,\delta)$-DP.
Moreover, if we set $\\{\gamma_{t}\\}=\gamma=\sqrt{T}$,
$T=\frac{n\epsilon}{M\lambda\sqrt{d\log(1/\delta)}}$ and
$\lambda=O((\frac{\sqrt{n\epsilon}M}{\kappa\sqrt[4]{d\log(1/\delta)}})^{\frac{2}{3}})$.
Then for any failure probability $\delta^{\prime}$, the output $\hat{w}$
satisfies the following with probability at least $1-\delta^{\prime}$

$\mathbb{E}[\mathcal{L}(\hat{w})]-L(\theta^{*})\leq
O\left(M^{\frac{4}{3}}\frac{\kappa^{\frac{2}{3}}(d\log(1/\delta))^{\frac{1}{6}}\sqrt{\log(1/\delta^{\prime})}}{(n\epsilon)^{\frac{1}{3}}}\right),$

where the expectation is taken over the randomness of noise, and the
probability is w.r.t. the dataset $D$.

###### Proof.

The proof of DP is just by the Generalizer Gaussian mechanism. For utility,
the proof is almost the same as in the proof for Theorem 10, while the only
difference is the noise. Similar to (27) we have the following result with
probability at least $1-2\exp(-\tau)$

$\mathcal{L}(\hat{w})-\mathcal{L}(\theta^{*})\leq O\left(\frac{\kappa
M^{2}\bar{\gamma}}{T}+\frac{M\sigma^{2}}{\lambda}+\frac{\sigma
M\sqrt{\tau}}{\sqrt{T}}+\frac{\sigma^{2}}{\bar{\gamma}}+\frac{M\sigma\sqrt{\tau}}{\sqrt{T}\bar{\gamma}}+\frac{dT^{2}\kappa(\beta
M+\lambda)^{2}\cdot\log(1/\delta)}{n^{2}\epsilon^{2}\bar{\gamma}}\right).$
(28)

Take $\bar{\gamma}=\sqrt{T}$ then we have

$\mathcal{L}(\hat{w})-\mathcal{L}(\theta^{*})\leq O\left(\frac{\kappa
M^{2}\sqrt{\tau}}{\sqrt{T}}+\frac{M^{2}}{\lambda}+\frac{dT^{3/2}\kappa\lambda^{2}\cdot\log(1/\delta)}{n^{2}\epsilon^{2}}\right).$

Take $T=\frac{n\epsilon}{M\lambda\sqrt{d\log(1/\delta)}}$ we have

$\mathcal{L}(\hat{w})-\mathcal{L}(\theta^{*})\leq O\left(\frac{\kappa
M\sqrt{\lambda}\sqrt[4]{d\log(1/\delta)}\sqrt{\tau}}{\sqrt{n\epsilon}}+\frac{M^{2}}{\lambda}\right).$

Take
$\lambda=O((\frac{\sqrt{n\epsilon}M}{\kappa\sqrt[4]{d\log(1/\delta)}})^{\frac{2}{3}})\geq\max\\{\beta,1\\}M$
we have w.p at least $1-\delta^{\prime}$

$\mathcal{L}(\hat{w})-\mathcal{L}(\theta^{*})\leq
O\left(M^{\frac{4}{3}}\frac{\kappa^{\frac{2}{3}}(d\log(1/\delta))^{\frac{1}{6}}\sqrt{\log(1/\delta^{\prime})}}{(n\epsilon)^{\frac{1}{3}}}\right).$

∎

###  B.7 Proof of Theorem 12

[23] study DP-SCO with heavy-tailed data in Euclidean space and propose an
$(\epsilon,\delta)$-DP algorithm for any $0<\epsilon,\delta<1$ that achieves
an expected excess population risk of
$O(M\frac{d}{\sqrt{n}}+\frac{\sqrt{M}d^{\frac{5}{4}}}{\sqrt{n\epsilon}})$,
where $M$ is the $\ell_{2}$-norm diameter of the constraint set $\mathcal{C}$,
under the following assumptions

######  Assumption 5.

1) The loss function $\ell(w,x)$ is non-negative, differentiable and convex
for all $w\in\mathcal{C}$. 2) The loss function is $\beta$-smooth. 3) The
gradient of $\mathcal{L}(w)$ at the optimum is zero. 4) There is a constant
$\sigma$ such that for all $j\in[d]$ and $w\in\mathcal{C}$ we have
$\mathbb{E}[\langle\nabla\ell(w,x)-\nabla\mathcal{L}(w),e_{j}\rangle^{2}]\leq\sigma^{2}$,
where $e_{j}$ is the $j$-th standard basis vector. 5) For any
$w\in\mathcal{C}$, the distribution of the gradient has bounded mean, i.e.,
$\|\nabla\mathcal{L}(w)\|_{2}\leq R$.

For $\ell_{p}^{d}$ space, we know that $L$-Lipschitz w.r.t $\|\cdot\|$ implies
$L$-Lipschitz w.r.t $\|\cdot\|_{2}$. Moreover,
$\mathbb{E}[||\nabla\ell(w,x)-\nabla\mathcal{L}(w)||_{*}^{2}]\leq\sigma^{2}$
implies
$\mathbb{E}[||\nabla\ell(w,x)-\nabla\mathcal{L}(w)||_{2}^{2}]\leq\sigma^{2}$
which indicates condition 4) in Assumption 5. For the diameter, it has the
diameter of $d^{\frac{1}{2}-\frac{1}{p}}M$ w.r.t $\|\cdot\|_{2}$. Thus we have
the following result.

######  Theorem 17.

For the $\ell_{p}^{d}$ space with $2\leq p\leq\infty$, suppose Assumption 4
holds. Then the Algorithm 1 in [23] is $(\epsilon,\delta)$-DP for all
$0<\epsilon,\delta<1$. Moreover, suppose the loss function is non-negative,
there exists $R=O(1)$ such that $\|\nabla\mathcal{L}(w)\|_{*}\leq R$ for all
$w\in\mathcal{C}$ and 3) in Assumption 5 holds. then the output satisfies

$\mathbb{E}[\mathcal{L}(w)]-\mathcal{L}(\theta^{*})\leq
O\left(\frac{d^{\frac{3}{2}-\frac{1}{p}}}{\sqrt{n}}+\frac{d^{\frac{3}{2}-\frac{1}{2p}}}{\sqrt{n\epsilon}}\right).$
(29)

∎


# Detection for Flat-Fading Channels based on Expectation Propagation

Elisa Conti, Amina Piemontese, Giulio Colavolpe, and Armando Vannucci The
authors are with the Department of Engineering and Architecture, University of
Parma, 43124 Parma, Italy (e-mails: {elisa.conti, amina.piemontese,
giulio.colavolpe, armando.vannucci}@unipr.it).

###### Abstract

This paper aims at tackling the problem of signal detection in flat-fading
channels. In this context, receivers based on the expectation propagation
framework appear to be very promising although presenting some critical
issues. We develop a new algorithm based on this framework where, unlike
previous works, convergence is achieved after a single forward-backward pass,
without additional inner detector iterations. The proposed message scheduling,
together with novel adjustments of the approximating distributions’
parameters, allows to obtain significant performance advantages with respect
to the state-of-the-art solution. Simulation results show the applicability of
this algorithm when sparser pilot configurations have to be adopted and a
considerable gain compared to the current available strategies.

###### Index Terms:

Expectation propagation, factor graphs, signal detection, fading channels,
parameter estimation, single-input single-output (SISO).

##  I Introduction

Detectors for flat Rayleigh fading channels have been widely studied in the
literature (e.g., see [1, 2, 3] and references therein). This kind of channels
is commonly adopted because it properly represents real scenarios, such as
narrowband communications over multipath wireless channels. The framework
based on factor graphs (FGs) and the sum-product algorithm (SPA) provides a
powerful tool for the accomplishment of both detection and decoding [4].
However, when this framework is applied jointly to detection and channel
estimation, discrete and continuous random variables appear in the FG and the
SPA becomes impractical.

A common approach foresees the use of canonical distributions [4]. Focusing on
the FG methodology and on the choice of a predefined family which allows some
kind of parametrization, different approaches have been analyzed [3, 1]. A
canonical distribution that is commonly selected for the messages in this
scenario is the Gaussian one. Due to the presence of both discrete and
continuous random variables, the observation message is typically represented
by a multimodal distribution. One of the main differences among the various
methods lies in the technique followed for its approximation. In [1], the
authors proposed to set the mean to the empirical average of the Gaussian
mixture and fix the variance to that of the additive white Gaussian noise
(AWGN). In [3], the approximating probability density function (pdf) is the
Gaussian that minimizes the Kullback-Leibler (KL) divergence from the true
message [5] and whose variance depends on the fading sample (i.e., on a
variable node), which complicates the resulting algorithm. Such a dependence
is removed by substituting the fading sample with its estimate performed at
the previous iteration. In the separate detection and decoding scenario, this
approach coincides with the one of [1]. By reducing the multimodal
distributions to single Gaussian ones, the forward and backward recursions of
the algorithm coincide with those of a Kalman smoother [6]. Although this
approach, based on message projection, is effective in pilot-aided
transmission, it eventually fails when the spacing between pilot blocks is
longer than a maximum accepted distance, as we will demonstrate in the
numerical results.

In order to overcome this limitation, in this paper we focus on a different
kind of receiver, which is based on the expectation propagation (EP) framework
proposed in [7] and generalized in [8]. This message-passing algorithm is
based on approximating the distributions of the variables of interest, rather
than those of the exchanged messages. This way, a prior belief coming from the
rest of the FG drives the approximation of intractable messages. The
combination of the SPA rules with the approximation deriving from EP leads to
a hybrid framework which, in the literature, is commonly referred to as EP-
belief propagation (BP). This framework has been deeply studied and analyzed
in different scenarios, including, e.g., the transmission over phase noise
channels [9], the MIMO sparse codes multiple access detection [10] or joint
channel estimation and decoding for faster-than-Nyquist signaling [11] over
frequency selective fading channels. As it is inherent to the EP update
strategy, improper distributions, e.g., Gaussians with a negative variance,
could appear. While in [2] the presence of improper distributions is accepted,
and no additional operations are introduced, different approaches are followed
by other authors, consisting in either rejecting the message updates that
produce improper distributions [12] or, in the case of a multivariate normal
distribution, operating on the covariance matrix through a replacement of
negative eigenvalues [13].

In [2], the EP framework has been applied to a flat-fading channel, modeled as
a complex autoregressive moving-average process. The authors proposed an EP
smoothing approach where inner detector iterations are required to achieve
convergence and it is tested using differential encoding.

However, when using a more powerful encoder, hence working in a lower signal-
to-noise ratio (SNR) regime, the approach in [2] is not equally effective and
this is the focus of our investigation.

The main contributions of this letter can be summarized as follows.

  * •

We propose an implementation of EP-BP that overcomes its limitations and
reduces the algorithmic complexity.

  * •

An improved scheduling which allows to achieve convergence without the aid of
multiple inner detector iterations (iterative refinement [7]) is adopted.

  * •

A new technique for balancing the variance of the involved Gaussian
distributions is introduced.

  * •

The problem of improper distributions is addressed, the scenarios where they
can compromise the algorithm operations are identified and a novel strategy to
handle these events is presented.

Finally, we demonstrate through simulations that the proposed techniques
present a superior performance, in terms of bit-error rate (BER), with respect
to state-of-the-art EP implementations [2]. Moreover, the proposed solution
ensures a higher robustness against concentrated pilot distributions, laying
the foundation for overcoming the limitations related to the classical and
well-established Kalman filter.

The rest of the paper is organized as follows. Section II introduces the
system model. Section III is devoted to the description of the EP algorithm.
The proposed techniques are introduced in Section IV. Numerical results are
presented in Section V while Section VI collects some concluding remarks.

##  II System Model

We consider the transmission of a sequence of $K$ complex coded symbols
$\boldsymbol{c}=\left(c_{0},c_{1},\dots,c_{K-1}\right)$ over a flat fading
channel. The sequence $\boldsymbol{c}$ is the result of the binary encoding,
through the code $\mathcal{C}$ and modulation over an $M$-ary constellation
($\mathcal{X}=\left\\{x_{1},\dots,x_{M}\right\\}\subset\mathbb{C}$) of the
information bit sequence $\boldsymbol{b}$. Assuming linear modulation with
Nyquist shaping pulses, matched filtering, and channel variations slow enough
so that no intersymbol interference arises, the discrete-time received signal
can be expressed as

$r_{k}=g_{k}c_{k}+n_{k},\quad k=0,1,\dots,K-1$ (1)

where ${n_{k}}$ is the AWGN and ${g_{k}}$ is the fading process. The vector of
noise samples, $\boldsymbol{n}=\left(n_{0},n_{1},\dots,n_{K-1}\right)$,111The
following notation is used throughout the paper. Bold lower-case letters
($\boldsymbol{a}$) denote vectors, bold upper-case letters ($\boldsymbol{A}$)
denote matrices. The superscripts $(\cdot)^{T}$ and $(\cdot)^{H}$ denote the
transpose and the Hermitian of a matrix or vector, respectively. Given a
square matrix $\boldsymbol{B}$ with dimension $N_{B}\times N_{B}$, the generic
$i,j$ element is denoted by $\boldsymbol{B}(i,j)$, $i,j=1,\dots,N_{B}$. has
independent and identically distributed (i.i.d) complex circularly symmetric
components with $n_{k}\sim\mathcal{N}_{\mathbb{C}}\left(0,N_{0}\right)$. 222A
complex circularly symmetric Gaussian random variable $v$ with mean $\mu$ and
variance $2\sigma^{2}$ is denoted by
$v\sim\mathcal{N}_{\mathbb{C}}\left(\mu,2\sigma^{2}\right)$. The corresponding
probability density function (pdf), with argument $x$, is denoted as
$g_{\mathbb{C}}\left(x,\mu,2\sigma^{2}\right)$. We adopt a similar notation
for a Gaussian vector, whose pdf is written as
$g_{\mathbb{C}}\left(\boldsymbol{x},\boldsymbol{\mu},\boldsymbol{C}\right)$,
where the inverse of the covariance matrix,
$\boldsymbol{\Sigma}=\boldsymbol{C}^{-1}$, is called precision.

The fading channel follows Clarke’s isotropic scattering model and the fading
coefficients $\boldsymbol{g}=\left(g_{0},g_{1},\dots,g_{K-1}\right)$ form a
sequence of zero-mean complex Gaussian random variables with unit variance,
statistically independent of both $\boldsymbol{c}$ and $\boldsymbol{n}$, with
autocorrelation sequence

$R_{g}(m)=E\left\\{{g_{k+m}g^{*}_{k}}\right\\}=J_{0}(2\pi f_{D}Tm)$ (2)

where $J_{0}(\cdot)$ is the zero-th order Bessel function, $T$ is the symbol
period, and $f_{D}$ is the Doppler spread. At the receiver, the adopted model
for the fading is an autoregressive process of order $N$, AR($N$), so that
$\boldsymbol{g}$ is approximated by the process
$\boldsymbol{h}=\left(h_{0},h_{1},\dots,h_{K-1}\right)$, which obeys333We
assume that $h_{l}=0$ for $l<0$.

$h_{k}=\sum_{n=1}^{N}\rho_{n}h_{k-n}+\nu_{k}=\boldsymbol{\rho}\boldsymbol{h}^{T}_{k-1}+\nu_{k}\,,$
(3)

where $\nu_{k}\sim\mathcal{N}_{\mathbb{C}}\left(0,2\sigma^{2}_{\nu}\right)$,
$\boldsymbol{h}_{k-1}=\left(h_{k-1},\dots,h_{k-N}\right)$ and
$\boldsymbol{\rho}=\left(\rho_{1},\rho_{2},\dots,\rho_{N}\right)$ are the real
AR($N$) coefficients. Clearly, there is a channel mismatch between the AR($N$)
process $\left\\{h_{k}\right\\}$ and the one modelled by (2), hence it is
critical to properly set the values of the coefficients $\boldsymbol{\rho}$ in
(3) and of the random increment’s variance $2\sigma^{2}_{\nu}$. A classical
approach in the literature [1, 14] consists in obtaining the coefficients
$\boldsymbol{\rho}$ by solving the Yule-Walker (YW) equations
$R_{h}(m)=R_{g}(m)$, for $m=1,2,\dots,N$, and computing $\sigma^{2}_{\nu}$ as
a function of $\boldsymbol{\rho}$ through the constraint $R_{h}(0)=R_{g}(0)$.
However, in the present context, it turns out that this solution is not the
best choice to approximate the actual fading statistics [3]. Therefore, we
will select the value of $\sigma^{2}_{\nu}$ that optimizes the performance, as
detailed in Section V.

##  III The Expectation Propagation Framework

Solving the fading estimation problem through message-passing algorithms
requires the factorization of the joint posterior pdf of the information bits
and the fading coefficients, given the received samples. From the system model
described above, it can be expressed as

$\begin{split}&P(\boldsymbol{b},\boldsymbol{c},\boldsymbol{h}|\boldsymbol{r})\propto
p(\boldsymbol{r}|\boldsymbol{c},\boldsymbol{h})p(\boldsymbol{h})P(\boldsymbol{c}|\boldsymbol{b})P(\boldsymbol{b})\\\
&\\!\propto\\!I[\boldsymbol{c}=\mu_{\mathcal{C}}(\boldsymbol{b})]p(h_{0})\\!\\!\prod_{k=1}^{K-1}\\!\\!p_{\nu}(h_{k}\\!-\boldsymbol{\rho}\boldsymbol{h}^{T}_{k-1})\\!\\!\prod_{k=0}^{K-1}f_{k}(c_{k},h_{k}),\end{split}$
(4)

where the information bits are uniformly and identically distributed,
$\mu_{\mathcal{C}}$ is the function that maps binary information messages
$\boldsymbol{b}$ into the codewords $\boldsymbol{c}$, and $I[\cdot]$ is an
indicator function equal to $1$ when
$\boldsymbol{c}=\mu_{\mathcal{C}}(\boldsymbol{b})$ and to zero otherwise. The
factorizations of the pdfs $p(\boldsymbol{h})$ and
$p(\boldsymbol{r}|\boldsymbol{c},\boldsymbol{h})$ in (4) come, respectively,
from the assumption that the fading process follows the AR($N$) model (3) and
from the fact that the AWGN channel is memoryless, given $\boldsymbol{h}$.

In (4), we defined

$p_{\nu}(h_{k}-\boldsymbol{\rho}\boldsymbol{h}^{T}_{k-1})\triangleq
p(h_{k}|\boldsymbol{h}_{k-1})=g_{\mathbb{C}}(h_{k},\boldsymbol{\rho}\boldsymbol{h}^{T}_{k-1},2\sigma^{2}_{\nu})\\\
$ (5) $f_{k}(c_{k},h_{k})\triangleq
p(r_{k}|c_{k},h_{k})=g_{\mathbb{C}}(r_{k},h_{k}c_{k},N_{0})$ (6)

as the functions associated with the factor nodes in the FG representation of
the joint pdf (4), which is reported in Fig. 1.

Figure 1: FG corresponding to (4).

Due to the presence of continuous random variables in the FG, we follow the
commonly adopted technique based on canonical distributions [4]. In this way,
only the parameters characterizing a specific pdf within the family are
exchanged in the FG. From (3), the distributions that better approximate those
of the fading samples are the Gaussian ones. According to the SPA rules [4],
which we assume the reader to be familiar with, the observation message, i.e.,
the message from factor node $f_{k}$ to variable node $\boldsymbol{h}_{k}$ is

$p_{d}(h_{k})\propto\sum^{M}_{m=1}P_{d}(c_{k}=x_{m})f_{k}(c_{k}=x_{m},h_{k})\,,$
(7)

hence it is the mixture of $M$ Gaussian distributions and its exact
propagation would be infeasible. Depending on how the message in (7) is
treated, different algorithms can be derived.

A straightforward projection of the mixture message (7) onto a single Gaussian
yields the Kalman filtering approach of [1, 3]. Once the projection is
performed, all the subsequent messages that are sent along the Markov chain at
the bottom of Fig. 1 belong to the same approximating family because of the
closure of the Gaussian family, under both the multiplication and the
convolution operations [1, 4].

A different possibility consists in approximating the distributions of
variable nodes, rather than the messages, and this is the key idea on which
the EP framework is based. The approximation of $p_{d}(h_{k})$ in (7) is thus
performed in EP, under the influence of a temporary prior distribution which
is the product of the incoming messages to the variable node
$\boldsymbol{h}_{k}$ from the lower part of the FG. When, on the contrary, a
message projection is carried out as in the Kalman filtering approach, so that
no prior information is exploited, this is sometimes called transparent
propagation algorithm [15, 16].

EP is a message-passing algorithm whose approximations are generated by
projecting the variable distribution on a fully-factorized exponential family
$\mathcal{F}$, where the projection is performed by minimizing the Kullback-
Leibler (KL) divergence. These two features allow to compute the projection by
solving a set of moment matching equations [8]. In the present case of a
Gaussian approximating family, the distribution resulting from the projection
is thus the unique Gaussian whose mean, variance and scale factor equal those
of the distribution to be approximated.

The specific algorithmic steps of EP are described in the following. Starting
from (7), the message $p_{d}(h_{k})$ can be rewritten as a function of the
channel parameter $h_{k}$ as follows

$p_{d}(h_{k})=\sum_{m=1}^{M}P_{m}\,g_{\mathbb{C}}(h_{k},\eta_{m},2\sigma^{2}_{m})\,,$
(8)

where

$P_{m}=\frac{P_{d}(c_{k}=x_{m})}{\left|x_{m}\right|^{2}}\hskip
5.69054pt,\hskip 5.69054pt\eta_{m}=\frac{r_{k}}{x_{m}}\hskip 5.69054pt,\hskip
5.69054pt2\sigma^{2}_{m}=\frac{N_{0}}{\left|x_{m}\right|^{2}}\hskip
5.69054pt.$ (9)

The dependence of $P_{m}$ and $\eta_{m}$ on the time epoch $k$ is omitted for
simplicity of notation. The temporary prior distribution
$p^{(p)}(\boldsymbol{h}_{k})$ is defined as

$\begin{split}&p^{(p)}(\boldsymbol{h}_{k})\triangleq
p^{(f)}_{k|k-1}(\boldsymbol{h}_{k})p^{(b)}_{k|k+1}(\boldsymbol{h}_{k})\\\
&=g_{\mathbb{C}}\left(\boldsymbol{h}_{k},\boldsymbol{\eta}_{k|k-1},\boldsymbol{C}_{k|k-1}\right)g_{\mathbb{C}}\left(\boldsymbol{h}_{k},\boldsymbol{\eta}_{k|k+1},\boldsymbol{C}_{k|k+1}\right)\\\
&\propto
g_{\mathbb{C}}\left(\boldsymbol{h}_{k},\boldsymbol{\eta}_{p,k},\boldsymbol{C}_{p,k}\right)\,,\end{split}$
(10)

where the $\propto$ sign derives from the messages normalization, and, by
defining $\boldsymbol{\Sigma}_{f,k}\triangleq\boldsymbol{C}_{k|k-1}^{-1}$ and
$\boldsymbol{\Sigma}_{b,k}\triangleq\boldsymbol{C}_{k|k+1}^{-1}$ as the
precisions of the forward and backward Gaussian messages, respectively, and
$\boldsymbol{\Sigma}_{p,k}\triangleq\boldsymbol{C}_{p,k}^{-1}$, the following
Gaussian product rules can be applied

$\displaystyle\boldsymbol{\Sigma}_{p,k}$
$\displaystyle=\boldsymbol{\Sigma}_{b,k}+\boldsymbol{\Sigma}_{f,k}\,,$ (11)
$\displaystyle{\boldsymbol{\Sigma}_{p,k}}\boldsymbol{\eta}_{p,k}$
$\displaystyle=\boldsymbol{\Sigma}_{f,k}\boldsymbol{\eta}_{k|k-1}+\boldsymbol{\Sigma}_{b,k}\boldsymbol{\eta}_{k|k+1}\,.$
(12)

The marginal distribution of $\boldsymbol{h}_{k}$ is given by the product of
all the incoming messages to the corresponding variable node, according to the
SPA rules [4]

$\begin{split}\hat{p}(\boldsymbol{h}_{k})&=p_{d}(\boldsymbol{h}_{k})p^{(f)}_{k|k-1}(\boldsymbol{h}_{k})p^{(b)}_{k|k+1}(\boldsymbol{h}_{k})\\\
&=\sum_{m=1}^{M}P_{m}\,g_{\mathbb{C}}\left(h_{k},\eta_{m},2\sigma^{2}_{m}\right)g_{\mathbb{C}}\left(\boldsymbol{h}_{k},\boldsymbol{\eta}_{p,k},\boldsymbol{C}_{p,k}\right)\\\
&=\sum_{m=1}^{M}P_{m}w_{m}\,g_{\mathbb{C}}\left(\boldsymbol{h}_{k},\boldsymbol{\eta}_{MIX,m},\boldsymbol{C}_{MIX,m}\right)\,,\end{split}$
(13)

where $w_{m}$ in (13) is the symbol-dependent weight that multiplies the
Gaussian distribution corresponding to the $m$-th symbol in (7) and is defined
as

$w_{m}=g_{\mathbb{C}}\left(\eta_{m},\boldsymbol{a}^{T}\boldsymbol{\eta}_{p,k},2\sigma^{2}_{m}+\boldsymbol{a}^{T}\boldsymbol{C}_{p,k}\boldsymbol{a}\right)\,,$
(14)

where $\boldsymbol{a}=\left(1,0,\dots,0\right)^{T}$ is a column vector of
dimension $N$ while $\boldsymbol{C}_{MIX,m}$ and $\boldsymbol{\eta}_{MIX,m}$,
after some calculations, here omitted for brevity, result:

$\displaystyle\boldsymbol{\eta}_{MIX,m}=\boldsymbol{\eta}_{p,k}+\boldsymbol{J}_{m,k}\left(\eta_{m}-\boldsymbol{a}^{T}\boldsymbol{\eta_{p,k}}\right)\,,$
(15)
$\displaystyle\boldsymbol{C}_{MIX,m}=\boldsymbol{C}_{p,k}-\boldsymbol{J}_{m,k}\,\boldsymbol{a}^{T}\boldsymbol{C}_{p,k}\,,$
(16)
$\displaystyle\boldsymbol{J}_{m,k}=\boldsymbol{C}_{p,k}\,\boldsymbol{a}\left(2\sigma^{2}_{m}+\boldsymbol{a}^{T}\boldsymbol{C}_{p,k}\,\boldsymbol{a}\right)^{-1}\,.$
(17)

Notice that the variable distribution in (13) is still a combination of $M$
Gaussians but the weight of each one, $w_{m}$, depends on the distance between
the mean of that mode ($\eta_{m}$) and that of the temporary prior
distribution. The closer they are, the more the corresponding mode is
reinforced by the prior and vice versa, as per (14).

The approximation of the mixture (13) with a simple Gaussian
$p^{(EP)}(\boldsymbol{h}_{k})=\text{proj}_{KL}\left[\hat{p}(\boldsymbol{h}_{k})\right]$
is performed, as usual in EP, by matching the moments of the distributions.
The zero-th order moment444We denote by $\mathbb{E}_{p}\left[x\right]$ the
expectation of $x$ under the distribution $p(\cdot)$.
$\mathbb{E}_{p^{(EP)}}\left[1\right]=\mathbb{E}_{\hat{p}}\left[1\right]$ is
simply the mass of the distribution, hence it would be necessary to set the
mass of $p^{(EP)}(\boldsymbol{h}_{k})$ to $Z=\sum_{m=1}^{M}P_{m}w_{m}$ that
results from (13). However, the mass of the messages is not significant in
message-passing algorithms like SPA or EP, therefore we prefer to consider
$p^{(EP)}(\boldsymbol{h}_{k})$ with unit mass (hence a proper pdf) and to
replace $\hat{p}(\boldsymbol{h}_{k})$ in (13) with its normalized version

$\\!\\!\\!p_{h}(\boldsymbol{h}_{k})\\!=\\!\frac{\hat{p}(\boldsymbol{h}_{k})}{Z}\\!=\\!\\!\sum_{m=1}^{M}\\!W_{m}g_{\mathbb{C}}\left(\boldsymbol{h}_{k},\boldsymbol{\eta}_{MIX,m},\boldsymbol{C}_{MIX,m}\right)\,.$
(18)

where $W_{m}=Z^{-1}P_{m}w_{m}$ are defined as normalized weights. The matching
of the first order moment is thus performed by equating the mean of the two
distributions
$\mathbb{E}_{p^{(EP)}}\left[\boldsymbol{h}_{k}\right]=\mathbb{E}_{p_{h}}\left[\boldsymbol{h}_{k}\right]$.
It follows that the mean of $p^{(EP)}(\boldsymbol{h}_{k})$ is

$\boldsymbol{\eta}_{EP}=\sum_{m=1}^{M}W_{m}\,\boldsymbol{\eta}_{MIX,m}\,.$
(19)

Then, by matching the second order moments, its covariance matrix is obtained
as

$\boldsymbol{C}_{EP}\\!=\\!\\!\sum_{m=1}^{M}\\!W_{m}\\!\left(\boldsymbol{C}_{MIX,m}+\boldsymbol{\eta}_{MIX,m}\boldsymbol{\eta}^{H}_{MIX,m}\right)-\boldsymbol{\eta}_{EP}\boldsymbol{\eta}^{H}_{EP}\,.$
(20)

Now, it is possible to compute the updated message

$\tilde{p}_{d}(\boldsymbol{h}_{k})\propto
g_{\mathbb{C}}\left(\boldsymbol{h}_{k},\boldsymbol{\eta}_{d},\boldsymbol{C}_{d}\right)\,,$
(21)

which is transmitted from the factor node $f_{k}$ to the variable node
$\boldsymbol{h}_{k}$ in Fig. 1, dividing the approximating Gaussian
$p^{(EP)}(\boldsymbol{h}_{k})$ by the temporary prior distribution
$p^{(p)}(\boldsymbol{h}_{k})$ in (10). We obtain

$\displaystyle\boldsymbol{\Sigma}_{d}$
$\displaystyle=\boldsymbol{\Sigma}_{EP}-\boldsymbol{\Sigma}_{p,k}\,,$ (22)
$\displaystyle\boldsymbol{\Sigma}_{d}\boldsymbol{\eta}_{d}$
$\displaystyle=\boldsymbol{\Sigma}_{EP}\boldsymbol{\eta}_{EP}-\boldsymbol{\Sigma}_{p,k}\boldsymbol{\eta}_{p,k}\,.$
(23)

The approximated message is then propagated following the SPA rules [4].

##  IV Proposed Algorithm

In the considered scenario, the classical EP implementation [2] presents some
critical issues which do not lead to a satisfactory performance. Therefore, in
the following, we will present novel techniques which allow to overcome the
limitations related to the standard EP structure.

###  IV-A Damping and Boosting Factors

In order to face the EP instabilities, as discussed in the literature, a
number of solutions exists, including damping [13, 12]. A dumped update of a
pdf is a member of the approximating family characterized by a convex
combination of new parameters (i.e., updated according to the current EP
iteration) and old ones (i.e., from previous iteration). However, unlike its
common use, in the analyzed scenario the introduction of this correction
factor is studied specifically to provide a performance improvement through
the balancing of overly confident estimates without the involvement of old
updates. In fact, in our implementation, the damping factor, $\xi_{d}$, which
belongs to the interval $[0,1]$, is applied to the precision of the
approximating message $\tilde{p}_{d}({h}_{k})$ only. In this way, the
algorithm underestimates the updates reliability as commonly done in the
literature when dealing with sub-optimal algorithms. We propose to set this
factor to the maximum value of $P_{d}(c_{k})$, i.e., the probability mass
function of the symbols. The main advantage of this approach is that a lower
number of turbo iterations is required for achieving convergence. At the first
turbo iteration555We will refer to the iterations between detector and decoder
as turbo iterations. (which is the only one in the separate detection and
decoding case) this quantity is equal to $1/M$. On the contrary, from the
$2$nd iteration ahead this value becomes a clearer and clearer representation
of the decoder confidence about the transmitted symbols. It follows that the
closer this value is to $1$, the higher the weight given to the approximating
message after the KL projection.

Moreover, we propose an additional technique which allows to obtain a further
performance improvement: the use of a boosting factor for pilot symbols.
Similarly to damping, the precision of the distribution in (21), whose
expression is reported in (22), is multiplied by a factor that, this time, is
greater than $1$, so that its effect is opposite: it increases the confidence
of the observation message in correspondence with a pilot symbol, i.e., when
$p_{d}({h}_{k})$ is a unimodal Gaussian distribution and, therefore, no
approximation is performed. The selected value for the boosting factor is
$\xi_{b}=2$.

###  IV-B Negative Variances

Due to the subtraction in (22), it may occur that the algorithm generates
negative variances, so that the approximation of the observation message,
$\tilde{p}_{d}({h}_{k})$, is not necessarily a valid probability distribution
[2]. This happens when the precision of the prior distribution is larger than
the one of the marginal after the KL projection. Since the focus of EP is on
the marginal distribution of the channel parameters, the algorithm could
accept improper distributions as messages. However, when the scheduling
described in Section IV-C is adopted for the separate detection and decoding
case, the only message which can present such variances is the observation
one.

When turbo iterations are considered, improper distributions could instead
propagate giving rise to numerical problems. Hence, in this case, we propose
an alternative approach to deal with negative variances, where the algorithm
performs a message projection instead of a variable distribution projection.
This leads to a significant performance improvement with respect to the
propagation of improper distributions as shown in Section V. The proposed
method consists in approximating the message $p_{d}(h_{k})$ locally with a
Gaussian distribution whose mean and variance are obtained through moment
matching without the influence of the rest of the network.

The approximating message can be written as

$\\!\hat{p}_{d}(h_{k})=\text{proj}_{\text{KL}}\left[\sum_{m=1}^{M}\frac{P_{d}(c_{m})}{\left|c_{m}\right|^{2}}\,g_{\mathbb{C}}\left(h_{k},\frac{r_{k}}{c_{m}},\frac{2\sigma^{2}_{n}}{\left|c_{m}\right|^{2}}\right)\right]\,.$
(24)

With the constraint on the Gaussian approximating family, the new mean and
variance will be

$\eta_{d}=\sum_{m=1}^{M}\frac{P_{d}(c_{k}=x_{m})}{\left|x_{m}\right|^{2}}\cdot\frac{r_{k}}{x_{m}}\,,\\\
$ (25)
$2\sigma^{2}_{d}=\sum_{m=1}^{M}\frac{P_{d}(c_{k}=x_{m})}{\left|x_{m}\right|^{2}}\cdot\left(\frac{N_{0}}{\left|x_{m}\right|^{2}}+\left|\frac{r_{k}}{x_{m}}\right|^{2}\right)-\left|\eta_{d}\right|^{2}\,.$
(26)

Therefore, if, following the main procedure where the KL divergence
minimization is performed on the marginal distribution, negative variances
occur, the resulting observation message update is discarded and the new
technique is applied.

###  IV-C Scheduling

We assume that the detector subgraph operates first by exchanging horizonal
messages along the Markov chain in the lower part of the FG shown in Fig. 1.
This message passing, implemented inside the detector, could be iterated with
multiple forward-backward passes ($N_{d}$) before sending vertical messages to
the upper part of the FG where the decoder operates. An inner message passing
procedure could take place also there until convergence or until a maximum
number of decoder iterations ($N_{c}$) is reached. The exchange of
information, through vertical messages, between detector and decoder can be
iterated for $N_{t}$ times.

With reference to the generic $n_{iter,T}$ turbo iteration, the proposed
algorithm structure is outlined in Algorithm 1.

Algorithm 1 Proposed EP algorithm

1:Initialization:

2:Initialize the messages $p^{(f)}_{0|-1}(\boldsymbol{h}_{0})$ and
$p^{(b)}_{K-1|K}(\boldsymbol{h}_{K-1})$ by setting
$\boldsymbol{\Sigma}_{f,0}$, $\boldsymbol{\Sigma}_{b,K-1}$ to $\mathbb{I}_{N}$
and $\boldsymbol{\eta}_{0|-1}$, $\boldsymbol{\eta}_{K-1|K}$ to
$\boldsymbol{0}$ according to the moments of the generated fading samples

3:for $k=0:K-2$ do

4: Compute the temporary prior distribution, i.e., update its parameters
$\boldsymbol{\eta}_{p,k}$ and $\boldsymbol{\Sigma}_{p,k}$ via (11)-(12). The
backward message $p^{(b)}_{k|k+1}(\boldsymbol{h}_{k})$ in (10) is set to the
one computed at the ($n_{iter,T}-1$)-th turbo iteration. If $n_{iter,T}=1$, it
is set to $1$ for any value of $k$

5: Compute the parameters of the exact marginal distribution (i.e., for each
$m$, compute the terms in (9), (14)-(17))

6: Perform the moment matching, using (19)-(20)

7: Compute $\tilde{p}_{d}(\boldsymbol{h}_{k})$ via (22)-(23) and perform
either damping or boosting on $\boldsymbol{\Sigma}_{d}$ depending on $k$

8: if $(\boldsymbol{\Sigma}_{d}(1,1)<0)\And(n_{iter,T}>1)$ then

9: Perform the observation message projection via (25)-(26)

10: end if

11: Compute the forward messages $p^{(f)}_{k|k}(\boldsymbol{h}_{k})$ and
$p^{(f)}_{k+1|k}(\boldsymbol{h}_{k+1})$ through the SPA rules [4] (refer to
Fig. 1)

12:end for

13:for $k=K-1:-1:1$ do

14: Compute $p^{(p)}(\boldsymbol{h}_{k})$ via (10)-(12) where
$p^{(f)}_{k|k-1}(\boldsymbol{h}_{k})$ is set to the one computed at the
($n_{iter,T}-1$)-th turbo iteration. If $n_{iter,T}=1$, it is set to $1$ for
all $k$s

15: Repeat the steps 5-10

16: Compute the backward messages $p^{(b)}_{k|k}(\boldsymbol{h}_{k})$ and
$p^{(b)}_{k-1|k}(\boldsymbol{h}_{k-1})$ through the SPA rules [4]

17:end for

18:Multiply the forward message obtained through 3-12 with the backward one
computed in 13-17 for each $k$

19:Compute for each $k$ and each $c_{k}\in\mathcal{X}$ the probabilities
$P_{u}(c_{k})$ to be sent to the decoder as $P_{u}(c_{k})\propto
g_{\mathbb{C}}\left(r_{k},c_{k}\boldsymbol{a}^{H}\boldsymbol{\eta}_{u,k},N_{0}+|c_{k}|^{2}\boldsymbol{a}^{H}\boldsymbol{C}_{u,k}\boldsymbol{a}\right)$,
where $\boldsymbol{\eta}_{u,k}$ and $\boldsymbol{C}_{u,k}$ derive from 18

Notice that the forward and backward recursions lead to two different
approximations for ${p}_{d}(h_{k})$. The main advantage of the adopted
scheduling is that a satisfying performance is already achieved after a single
forward-backward sweep inside the detector.

##  V Simulation Results

The performance of the proposed algorithm was analyzed in terms of BER versus
$E_{b}/N_{0}$, $E_{b}$ being the received signal energy per information bit
and $N_{0}$ the one-sided noise power spectral density. We considered a
$(3,6)$-regular rate-$1/2$ LDPC code with length $4000$ and QPSK modulation
($M=4$). In all simulations, in order to allow the algorithms convergence,
pilot symbols were inserted in the transmitted codeword. 666We use the
notation $\\{$number of pilots in a block$\\}/\\{$blocks distance$\\}$ to
describe a pilot distribution. The presence of pilots involves a slight
decrease of the effective information rate, resulting in an increase in the
required SNR. This increase has been artificially introduced in the curve
labeled “All Pilots” , inserted as ideal bound, for the sake of comparison.
Hence, the gap between the “All Pilots” and the other curves is uniquely due
to the need for fading estimation/compensation, and not to the decrease due to
pilot symbols. The curves labeled “Kalman filter” in Figs. 2 and 4 refer to
the algorithm proposed in [3].

Figure 2: Performance for the $1/20$ pilots distribution with
$f_{D}T=10^{-2}$.

The algorithms were tested for both the AR($1$) and AR($2$) fading models. The
values of the parameter $\boldsymbol{\rho}$ were set according to the YW
equations, whereas the increment variance $\sigma^{2}_{\nu}$ was set to the
one minimizing the BER, found through computer simulations.

In Fig. 2 we compared the performance of the proposed algorithm with both the
Kalman filter and the “batch-EP” of [2], labeled with “classical EP”, when the
normalized Doppler bandwidth of the fading process $f_{D}T$ was set to
$10^{-2}$ and one pilot symbol was inserted in every block of $20$ transmitted
symbols. Firstly, we analyzed the practical case of separate detection and
decoding.777In order to indicate the number of inner detector, inner decoder
and turbo iterations, the notation $N_{d}$-$N_{c}$-$N_{t}$ is used. Figure 2
shows that the proposed EP algorithm can achieve a performance gain of
approximately $0.9\,$dB with respect to the Kalman filter for both the AR($1$)
and the AR($2$) models. Moreover, a significant improvement is observed
comparing the proposed technique with the classical EP of [2]. Regarding the
latter algorithm, the authors proposed an EP smoothing approach where an
iterative refinement of the observation approximations is foreseen.
Consistently with that used in [2], a maximum number ($N_{d}$) of $5$ inner
detector iterations was allowed, and it resulted to be sufficient for
achieving convergence. For this algorithm, the curve shown in Fig. 2 is
obtained using the AR($1$) model since the AR($2$) model did not yield further
improvement. The complexity of the considered algorithms was evaluated with
reference to one inner detector iteration for a single time epoch $k$ and it
is reported in Fig. 3. They were compared considering the total number of
required operations, including two-terms multiplications and additions. The
distinction among the various kind of operations is, instead, shown in Table I
for $M=4$ where the number of Lookup Table (LUT) accesses for the computation
of nonlinear terms is also specified. It can be noticed that the additional
operations involved by the EP framework lead inevitably to a complexity
increase which, relying on our implementation, is still limited and reasonably
comparable to the one of the Kalman filter. On the other hand, a considerable
complexity increase is registered for the classical strategy [2].

Figure 3: Detection algorithms complexity comparison for $N_{d}=1$. | Two-terms additions | Two-terms multiplications | LUT accesses
---|---|---|---
|  $\boldsymbol{N=1}$ | $\boldsymbol{N=2}$ |  ${\boldsymbol{N=1}}$ | ${\boldsymbol{N=2}}$ | $-$
Kalman | ${\boldsymbol{24}}$ | ${\boldsymbol{97}}$ | ${\boldsymbol{53}}$ | ${\boldsymbol{177}}$ | $\boldsymbol{6}$
EP proposed | ${\boldsymbol{86}}$ | ${\boldsymbol{261}}$ | ${\boldsymbol{148}}$ | ${\boldsymbol{452}}$ | $\boldsymbol{12}$
EP classical | ${\boldsymbol{217}}$ | ${\boldsymbol{602}}$ | ${\boldsymbol{478}}$ | ${\boldsymbol{1138}}$ | $\boldsymbol{8}$
TABLE I:  Complexity comparison for $M=4$ and $N_{t}=1$.

In Fig. 2, the BER curves related to the joint iterative detection and
decoding, where a maximum of $200$ turbo iterations is allowed, are also
reported with reference to the AR($1$) fading model. It is possible to notice
that the curves related to the proposed EP algorithm and the Kalman filter
(optimized according to [3]) are approximately overlapped. In fact, when a
sufficient number of turbo iterations is performed, the decoder provides
symbol probabilities approaching the ones passed in correspondence with a
known symbol. Therefore, the observation message is composed by a single
Gaussian and the operations performed by the two detectors coincide.

Figure 4: Performance for different pilots distribution with $f_{D}T=5\cdot
10^{-3}$ and $\sigma^{2}_{\nu}=8\cdot 10^{-3}$.

In Fig. 4, different pilot distributions were considered in order to show the
higher robustness of EP with respect to the Kalman filter against pilot
symbols arrangements in sequences separated by longer blocks of code symbols.
In particular, the analyzed configuration are $1/20$, $2/40$, $4/80$, $8/160$,
so that the effective information rate is the same. The normalized Doppler
bandwidth was set to $5\cdot 10^{-3}$. The AR($1$) fading model was adopted
for the sake of simplicity, however, similar conclusions can be drawn for
AR($2$). It is possible to notice that both the proposed algorithm and the
Kalman filter are practically insensitive to a change in the pilots blocks
distance from $20$ to $40$ symbols. Differently, increasing the gap to $80$
symbols causes a loss of about $0.8\,$dB for the EP-based algorithm and about
$1.5\,$dB for the Kalman filter approach. Finally, considering a pilots
distance of $160$ code symbols, the EP detector is still able to provide
reasonable fading estimates and, with the aid of $6$ turbo iterations, it
loses approximately $2\,$dB from the case where the $1/20$ pilots distribution
is adopted. On the contrary, in the same conditions, the Kalman filter does
not work because of its lower robustness against longer pilots spacing.

The proposed solution in this last situation adopts the observation message
projection strategy presented in Section IV-B. In order to show the advantages
of this approach, we compared the algorithm performance with and without the
handling of improper distributions in Fig. 5. It can be noted that a
considerable gain can be achieved by blocking the improper distributions
propagation. This is due to the fact that the adopted pilots configuration
makes the algorithm operating conditions harder. Therefore, the estimate
obtained by the prior is not often in accordance with any of the observation
message modes, leading to a rather wide approximating marginal distribution.
Negative variances are thus typically generated after the observation message
update, and allowing their propagation and combination with those from
previous iterations leads to poor performance. On the other hand, by
considering more favourable conditions, the performance gap becomes less
relevant.

Figure 5: Performance comparison of the EP algorithm adopting the strategy
proposed in Section IV-B with the procedure based on the propagation of
improper distributions. The analyzed scenario is based on the concentrated
pilots configuration, in particular the $8/160$ one, and the normalized
Doppler bandwidth $f_{d}T=5\cdot 10^{-3}$. The selected $\sigma^{2}_{\nu}$
optimizing the performance is set to $8\cdot 10^{-3}$. The simulated
iterations are $1$-$200$-$7$.

##  VI Conclusions

In this paper, we proposed a new solution for signal detection over channels
with flat-fading through EP. The presented method allows to solve the
numerical instabilities related to the traditional procedure provided in the
literature. We showed through numerical results the advantages in terms of BER
that can be achieved in the separate detection and decoding case with respect
to the classical Kalman filter approach and we demonstrated the higher
robustness of this algorithm against more concentrated pilots distributions.

The advantages of this framework offer a variety of possibilities for new
research directions, such as its application to the multiple-input multiple-
output channel model and the study of more severe fading conditions, for
instance frequency selectivity. Moreover, we could extend our analysis to
other kinds of divergence measures for the involved approximations.

## References

  * [1] H. Niu, M. Shen, J. A. Ritcey, and H. Liu, “A factor graph approach to iterative channel estimation and LDPC decoding over fading channels,” _IEEE Trans. Wireless Commun._ , vol. 4, no. 4, pp. 1345–1350, Jul. 2005\.
  * [2] Y. Qi and T. P. Minka, “Window-based expectation propagation for adaptive signal detection in flat-fading channels,” _IEEE Trans. Wireless Commun._ , vol. 6, no. 1, pp. 348–355, 2007.
  * [3] A. Barbieri, A. Piemontese, and G. Colavolpe, “On the ARMA approximation for fading channels described by the Clarke model with applications to Kalman-based receivers,” _IEEE Trans. Wireless Commun._ , vol. 8, no. 2, pp. 535–540, 2009.
  * [4] F. R. Kschischang, B. J. Frey, and H.-A. Loeliger, “Factor graphs and the sum-product algorithm,” _IEEE Trans. Inform. Theory_ , vol. 47, pp. 498–519, Feb. 2001.
  * [5] T. M. Cover and J. A. Thomas, _Elements of Information Theory_. New York: John Wiley & Sons, Inc., 1991.
  * [6] M. I. Jordan and C. Bishop, “An introduction to graphical models,” 2001.
  * [7] T. P. Minka, “Expectation propagation for approximate Bayesian inference,” in _Proc. UAI_ , Aug. 2001.
  * [8] ——, “Divergence measures and message passing,” Microsoft Research Cambridge, Tech. Rep., 2005, mSR-TR-2005-173 (2005).
  * [9] G. Colavolpe, E. Conti, A. Piemontese, and A. Vannucci, “The difficult road of expectation propagation towards phase noise detection,” to be presented at IEEE Intern. Conf. on Commun., Rome (Italy), May 2023.
  * [10] W. Yuan, N. Wu, Q. Guo, Y. Li, C. Xing, and J. Kuang, “Iterative receivers for downlink MIMO-SCMA: Message passing and distributed cooperative detection,” _IEEE Transactions on Wireless Communications_ , vol. 17, no. 5, pp. 3444–3458, 2018.
  * [11] N. Wu, W. Yuan, Q. Guo, and J. Kuang, “A hybrid BP-EP-VMP approach to joint channel estimation and decoding for ftn signaling over frequency selective fading channels,” _IEEE Access_ , vol. 5, pp. 6849–6858, 2017.
  * [12] M. Seeger, “Bayesian Gaussian process models: PAC-bayesian generalisation error bounds and sparse approximations,” _PhD thesis, University of Edinburgh_ , 2003.
  * [13] A. Gelman, A. Vehtari, P. Jylänki, T. Sivula, D. Tran, S. Sahai, P. Blomstedt, J. P. Cunningham, D. Schiminovich, and C. Robert, “Expectation propagation as a way of life: A framework for bayesian inference on partitioned data,” _arXiv preprint arXiv:1412.4869_ , 2014.
  * [14] J. Heo and K. M. Chugg, “Adaptive iterative detection for turbo codes on flat fading channels,” in _WCNC_ , Chicago, IL, Sep. 2000.
  * [15] A. Vannucci, G. Colavolpe, R. Pecori, and L. Veltri, “Estimation of a Gaussian source with memory in bursty impulsive noise,” in _2019 IEEE ISPLC_ , 2019, pp. 1–6.
  * [16] A. Vannucci, G. Colavolpe, and L. Veltri, “Estimation of correlated Gaussian samples in impulsive noise,” _IEEE Commun. Letters_ , vol. 24, no. 1, pp. 103–107, 2020.