Multi-User Detection Based on Expectation Propagation for the Non-Coherent SIMO Multiple Access Channel

TL;DR

This paper introduces a novel expectation propagation-based multi-user detector for non-coherent SIMO multiple access channels, achieving improved performance over existing methods in symbol error rate and sum-rate, with polynomial complexity.

Contribution

It develops a new EP-based detector for non-coherent SIMO channels that outperforms state-of-the-art methods and simplifies complexity with two variants.

Findings

Outperforms projection-based interference mitigation methods.

Achieves lower symbol error rates than conventional coherent detectors.

Provides high achievable sum-rates and improved BER with turbo decoding.

Abstract

We consider the non-coherent single-input multiple-output (SIMO) multiple access channel with general signaling under spatially correlated Rayleigh block fading. We propose a novel soft-output multi-user detector that computes an approximate marginal posterior of each transmitted signal using only the knowledge about the channel distribution. Our detector is based on expectation propagation (EP) approximate inference and has polynomial complexity in the number of users, the number of receive antennas, and channel coherence time. We also propose two simplifications of this detector with reduced complexity. With Grassmannian signaling, the proposed detectors outperform a state-of-the-art non-coherent detector with projection-based interference mitigation. With pilot-assisted signaling, the EP detector outperforms, in terms of symbol error rate, some conventional coherent pilot-based detectors, including a sphere decoder and a joint channel estimation-data detection scheme. Our EP-based detectors produce accurate approximates of the true posterior leading to high achievable sum-rates. The gains of these detectors are further observed in terms of the bit error rate when using their soft outputs for a turbo channel decoder.

Equations

$$\displaystyle{\mathcal{N}}({\boldsymbol{x}};{\boldsymbol{\mu}},\boldsymbol{\Sigma}):=\frac{\exp\big{(}-({\boldsymbol{x}}-{\boldsymbol{\mu}})^{\scriptscriptstyle\mathsf{H}}\boldsymbol{\Sigma}^{-1}({\boldsymbol{x}}-{\boldsymbol{\mu}})\big{)}}{\pi^{n}{\rm det}\left(\right)(\boldsymbol{\Sigma})},~{}{\boldsymbol{x}}\in\mathbb{C}^{n}$$

$$\mathsf{Y}$$

$$p_{{\boldsymbol{\mathsf{Y}}}|{\boldsymbol{\mathsf{S}}}}({\boldsymbol{Y}}|{\boldsymbol{S}})=\\ \frac{\exp\big{(}-{\rm vec}({\boldsymbol{Y}}^{\scriptscriptstyle\mathsf{T}})^{\scriptscriptstyle\mathsf{H}}\big{(}\sigma^{2}{\boldsymbol{I}}_{NT}+\sum_{k=1}^{K}{\boldsymbol{s}}_{k}{\boldsymbol{s}}_{k}^{\scriptscriptstyle\mathsf{H}}\otimes\boldsymbol{\Xi}_{k}\big{)}^{-1}{\rm vec}({\boldsymbol{Y}}^{\scriptscriptstyle\mathsf{T}})\big{)}}{\pi^{NT}{\rm det}\left(\right)(\sigma^{2}{\boldsymbol{I}}_{NT}+\sum_{k=1}^{K}{\boldsymbol{s}}_{k}{\boldsymbol{s}}_{k}^{\scriptscriptstyle\mathsf{H}}\otimes\boldsymbol{\Xi}_{k})}.$$

$$\hat{\mathbf{S}}$$

$${\mathrm{LLR}}_{k,j}(\mathbf{Y})$$

$$=\log \frac{{\sum }_{\mathbf{\alpha }\in {\mathcal{S}}_{k,j}^{(1)}}{p}_{\mathsf{Y}|{\mathsf{s}}_k}(\mathbf{Y}|\mathbf{\alpha })}{{\sum }_{\mathbf{\beta }\in {\mathcal{S}}_{k,j}^{(0)}}{p}_{\mathsf{Y}|{\mathsf{s}}_k}(\mathbf{Y}|\mathbf{\beta })}$$

$$=\log \frac{{\sum }_{\mathbf{\alpha }\in {\mathcal{S}}_{k,j}^{(1)}}{p}_{{\mathsf{s}}_k|\mathsf{Y}}(\mathbf{\alpha }|\mathbf{Y})}{{\sum }_{\mathbf{\beta }\in {\mathcal{S}}_{k,j}^{(0)}}{p}_{{\mathsf{s}}_k|\mathsf{Y}}(\mathbf{\beta }|\mathbf{Y})}$$

$$p_{\mathsf{S}|\mathsf{Y}}(\mathbf{S}|\mathbf{Y})=\frac{p_{\mathsf{Y}|\mathsf{S}}(\mathbf{Y}|\mathbf{S})p_{\mathsf{S}}(\mathbf{S})}{p_{\mathsf{Y}}(\mathbf{Y})}\propto p_{\mathsf{Y}|\mathsf{S}}(\mathbf{Y}|\mathbf{S})p_{\mathsf{S}}(\mathbf{S}).$$

$$\mathrm{Pr}(\mathsf{S}=[{\mathbf{s}}_1,\dots ,{\mathbf{s}}_K])=\prod _{k=1}^K\frac{1}{|{\mathcal{S}}_k|}\text{mathbbm}1\{{\mathbf{s}}_k\in {\mathcal{S}}_k\}.$$

$$p_{{\mathsf{s}}_k|\mathsf{Y}}({\mathbf{s}}_k|\mathbf{Y})=\sum _{{\mathbf{s}}_l\in {\mathcal{S}}_l,\forall l\ne k}{p}_{\mathsf{S}|\mathsf{Y}}([{\mathbf{s}}_1,\dots ,{\mathbf{s}}_K]|\mathbf{Y}),k\in [K].$$

$$p_{\mathsf{S}|\mathsf{Y}}([{\mathbf{s}}_1,\dots ,{\mathbf{s}}_K]|\mathbf{Y})$$

$$=\prod _{k=1}^K{\hat{p}}_{{\mathsf{s}}_k|\mathsf{Y}}({\mathbf{s}}_k|\mathbf{Y}).$$

$${\mathbf{C}}_{k1}$$

$${\mathbf{\mu }}_{k1}$$

$${\mathbf{C}}_{k1}(t)$$

$${\mathbf{\mu }}_{k1}(t)$$

$$\displaystyle\quad\times\bigg{(}\boldsymbol{\Sigma}_{k}(t)-\sum_{i=1}^{|{\mathcal{S}}_{k}|}\pi_{k1}^{(i)}(t)\boldsymbol{\Sigma}_{ki}(t)\bigg{)}{\boldsymbol{C}}_{k0}^{-1}(t\!-\!1){\boldsymbol{\mu}}_{k0}(t\!-\!1)$$

$$+(1-\eta ){\mathbf{\mu }}_{k1}(t-1),$$

$${\mathbf{C}}_{l0}(t)$$

$$\forall l\ne k,$$

$${\mathbf{\mu }}_{l0}(t)$$

$${\mathbf{s}}_k^{(i)}=\frac{{\mathbf{U}}_k{\mathbf{d}}^{(i)}}{\Vert {\mathbf{U}}_k{\mathbf{d}}^{(i)}\Vert },i\in [2^B].$$

$$\underset{{\mathbf{d}}^{(i)}\in G({\mathbb{C}}^{T-K+1},1),i=1,\dots ,2^B}{\max }\underset{1\le i<j\le 2^B}{\min }d({\mathbf{d}}^{(i)},{\mathbf{d}}^{(j)}),$$

$${\Delta }_t=\frac{1}{K}\sum _{k=1}^K\mathbb{E}{\left[\right]}_{\mathsf{Y}}[\mathcal{T}\mathcal{V}({\hat{p}}_{{\mathsf{s}}_k|\mathsf{Y}}^{(t)},p_{{\mathsf{s}}_k|\mathsf{Y}})]$$

$$\displaystyle D\big{(}q_{\alpha}({\boldsymbol{x}})\big{\|}\underline{p}({\boldsymbol{x}})\big{)}=\int q_{\alpha}({\boldsymbol{x}})\ln\frac{q_{\alpha}({\boldsymbol{x}})}{\prod_{\beta}\underline{p}_{\beta}({\boldsymbol{x}}_{\beta})}\mathop{}\mathrm{d}{\boldsymbol{x}}$$

$$=\sum _{\beta }\int q_{\alpha }(\mathbf{x})\ln \frac{1}{{\underline{p}}_{\beta }({\mathbf{x}}_{\beta })}\mathrm{d}\mathbf{x}+c_0$$

$$=\sum _{\beta \in {\mathfrak{N}}_{\alpha }}\int q_{\alpha }(\mathbf{x})\ln \frac{1}{{\underline{p}}_{\beta }({\mathbf{x}}_{\beta })}\mathrm{d}\mathbf{x}+\sum _{\beta \notin {\mathfrak{N}}_{\alpha }}\int q_{\alpha }(\mathbf{x})\ln \frac{1}{{\underline{p}}_{\beta }({\mathbf{x}}_{\beta })}\mathrm{d}\mathbf{x}$$

$$+c_0$$

$$=\sum _{\beta \in {\mathfrak{N}}_{\alpha }}\int q_{\alpha }(\mathbf{x})\ln \frac{1}{{\underline{p}}_{\beta }({\mathbf{x}}_{\beta })}\mathrm{d}\mathbf{x}$$

$$+\sum _{\beta \notin {\mathfrak{N}}_{\alpha }}\int {\hat{p}}_{\beta }({\mathbf{x}}_{\beta })\ln \frac{1}{{\underline{p}}_{\beta }({\mathbf{x}}_{\beta })}\mathrm{d}{\mathbf{x}}_{\beta }+c_0$$

$$\displaystyle=-\!\sum_{\beta\in\mathfrak{N}_{\alpha}}\!\int\!q_{\alpha}({\boldsymbol{x}})\left[\underline{{\boldsymbol{\gamma}}}_{\beta}^{\scriptscriptstyle\mathsf{T}}{\boldsymbol{\phi}}({\boldsymbol{x}}_{\beta})\!-\!A_{\beta}(\underline{{\boldsymbol{\gamma}}}_{\beta})\right]\mathop{}\mathrm{d}{\boldsymbol{x}}+\!\sum_{\beta\notin\mathfrak{N}_{\alpha}}\!D\big{(}\hat{p}_{\beta}\big{\|}\underline{p}_{\beta}\big{)}$$

$$+c_0$$

Full text

I Introduction

In wireless communications, multi-antenna based multiple-input multiple-output (MIMO) technology is capable of improving significantly both the system spectral efficiency and reliability due to its multiplexing and diversity gains [2, 3]. MIMO is at the heart of current cellular systems, and large-scale (massive) MIMO [4] is considered as one of the fundamental technologies for the fifth-generation (5G) wireless communications [5]. In practical MIMO systems, the transmitted symbols are normally drawn from a finite discrete constellation to reduce complexity. Due to propagation effects, the symbols sent from different transmit antennas interfere, and the receiver observes a linear superposition of these symbols corrupted by noise. The task of the receiver is to detect these symbols (or rather the underlying bits) based on the received signal and the available knowledge about the channel.

If the instantaneous value of the channel matrix is treated as known, the detection is said to be coherent and has been investigated extensively in the literature [6]. In this case, the data symbols are normally taken from a scalar constellation such as the quadrature amplitude modulation (QAM). Since the optimal maximum-likelihood (ML) coherent detection problem is known to be non-deterministic polynomial-time hard (NP-hard) [7], many sub-optimal coherent MIMO detection algorithms have been proposed. These range from linear schemes, such as the zero forcing (ZF) and minimum mean square error (MMSE) detectors, to non-linear schemes based on, for example, interference cancellation, tree search, and lattice reduction [6].

If only statistical information about the channel is available, the detection problem is said to be non-coherent. In the block fading channel where the channel matrix remains constant for each length-$T$ coherence block and varies between blocks, the receiver can estimate (normally imperfectly) the channel based on the transmitted pilot symbols, then perform coherent detection based on the channel estimate. Channel estimation and data detection can also be done iteratively [8, 9], or jointly based on tree search [10, 11]. These schemes requires pilot transmission for an initial channel estimate or to guarantee the identifiability of the data symbols. Another approach not involving pilot transmission is unitary space time modulation, in which the matrix of symbols in the space-time domain is orthonormal and isotropically distributed [12]. There, information is carried by the signal matrix subspace position, which is invariant to multiplication by the channel matrix. Therefore, a constellation over matrix-valued symbols can be designed as a collection of subspaces in $\mathbb{C}^{T}$. Such constellations belong to the Grassmann manifold $G(\mathbb{C}^{T},K)$, which is the space of $K$-dimensional subspaces in $\mathbb{C}^{T}$, where $K$ is the number of transmit antennas. For the independent and identically distributed (i.i.d.) Rayleigh block fading channel, when the signal-to-noise-ratio (SNR) is large, Grassmannian signaling was shown to achieve a rate within a vanishing gap from the capacity if $T\geq N+\min\{K,N\}$ [13], and within a constant gap if $2K\leq T\leq N+K$ [14], where $N$ is the number of receive antennas. Like with coherent detection, the optimal ML non-coherent detection problem under Grassmannian signaling is NP-hard. Thus, low-complexity sub-optimal detectors have been proposed for constellations with additional structure, e.g., [15, 16, 17].

In this paper, we focus on the non-coherent detection problem in the Rayleigh flat and block fading single-input multiple-output (SIMO) multiple-access channel (MAC) with coherence time $T$. There, the communication signals are independently transmitted from $K$ single-antenna users. If the users could cooperate, the high-SNR optimal joint signaling scheme would be a Grassmannian signaling on $G(\mathbb{C}^{T},K)$ [13]. However, we assume uncoordinated users, for which the optimal non-coherent transmission scheme is not known, although some approximate optimality design criteria have been proposed in [18]. In this work, we design the detector without assuming any specific structure of the signal transmitted over a coherence block. We consider the case where the receiver is interested not only in the hard detection of the symbols but also in their posterior marginal probability mass functions (PMFs). This “soft” information is needed, for example, when computing the bit-wise log-likelihood ratios (LLRs) required for soft-input soft-output channel decoding. Computing an exact marginal PMF would require enumerating all possible combinations of other-user signals, which is infeasible with many users, many antennas, or large constellations. Thus, we seek sub-optimal schemes with practical complexity.

In contrast to probabilistic coherent MIMO detection, for which many schemes have been proposed (e.g., [19, 20, 21]), the probabilistic non-coherent MIMO detection under general signaling, and Grassmannian signaling in particular, has not been well investigated. The detection scheme proposed in [22] is sub-optimal and compatible only with the specific constellation structure considered therein. The list-based soft demapper in [23] reduces the number of terms considered in posterior marginalization by including only those symbols at a certain distance from a reference point. However, it was designed for the single-user case only and has no obvious generalization to the MAC. The semi-blind approaches [8, 9, 10, 11] for the MIMO point-to-point channel can be extended to the MAC. However, these schemes are restricted to transmitted signals with pilots.

In this work, we propose message-passing algorithms for posterior marginal inference of non-coherent multi-user MIMO transmissions over spatially correlated Rayleigh block fading channels. Our algorithms are based on expectation propagation (EP) approximate inference [24, 25]. EP provides an iterative framework for approximating posterior beliefs by parametric distributions in the exponential family [26, Sec. 1.6]. Although there are many possible ways to apply EP to our non-coherent multi-user detection problem, we do so by choosing as variable nodes the indices of the transmitted symbols and the noiseless received signal from each user. The EP algorithm passes messages between the corresponding variable nodes and factor nodes on a bipartite factor graph. In doing so, the approximate posteriors of these variables are iteratively refined. We also address numerical implementation issues.

To measure the accuracy of the approximate posterior generated by the soft detectors, we compute the mismatched sum-rate of the system that uses the approximate posterior as the decoding metric. This mismatched sum-rate approaches the achievable rate of the system as the approximate posterior gets close to the true posterior. We also evaluate the symbol error rate when using the proposed schemes for hard detection, and the bit error rate when using these schemes for turbo equalization with a standard turbo code.

The contributions of this work are summarized as follows:

We propose soft and hard multi-user detectors for the non-coherent SIMO MAC using EP approximate inference, and methods to stabilize the EP updates. The proposed detectors work for general vector-valued transmitted symbols within each channel coherence block, i.e., it is general enough to include both the pilot-assisted and pilot-free signaling cases. 2. 2.

We propose two simplifications of the EP detector with reduced complexity. The first one, so-called EPAK, is based on approximating the EP messages with Kronecker products. The second one can be interpreted as soft MMSE estimation and successive interference approximation (SIA). 3. 3.

We analyze the complexity and numerically evaluate the convergence, running time, and performance of the proposed EP, EPAK, and MMSE-SIA detectors, the optimal ML detector, a genie-aided detector, the state-of-the-art detector from [22], and some conventional coherent pilot-based schemes. Our results suggest that the proposed detectors offer significantly improved mismatched sum-rate, symbol error rate, and coded bit error rate with respect to (w.r.t.) some existing sub-optimal schemes, while having lower complexity than the ML detector.

To the best of our knowledge, our proposed approach is the first message-passing scheme for non-coherent multi-user MIMO detection with general constellations.

The remainder of this paper is organized as follows. The system model is presented in Section II. A brief review of EP is presented in Section LABEL:sec:EP, and the EP approach to non-coherent detection is presented in Section LABEL:sec:EP_application. In Section LABEL:sec:MMSE-SIA_EPAK, two simplifications (MMSE-SIA and EPAK) of the EP detector are presented. Implementation aspects of EP, MMSE-SIA, and EPAK are discussed in Section LABEL:sec:implementation. Numerical results and conclusions are presented in Section VII and Section VIII, respectively. The mathematical preliminaries and proofs are provided in the appendices.

Notations: Random quantities are denoted with non-italic letters with sans-serif fonts, e.g., a scalar ${\mathsf{x}}$, a vector ${\boldsymbol{\mathsf{v}}}$, and a matrix ${\boldsymbol{\mathsf{M}}}$. Deterministic quantities are denoted with italic letters, e.g., a scalar $x$, a vector $\boldsymbol{v}$, and a matrix $\boldsymbol{M}$. The Euclidean norm is denoted by $\|{\boldsymbol{v}}\|$ and the Frobenius norm $\|{\boldsymbol{M}}\|_{F}$. The conjugate, transpose, conjugate transpose, trace, and vectorization of ${\boldsymbol{M}}$ are denoted by ${\boldsymbol{M}}^{*}$, ${\boldsymbol{M}}^{\scriptscriptstyle\mathsf{T}}$, ${\boldsymbol{M}}^{\scriptscriptstyle\mathsf{H}}$, ${\rm tr}\left\{\right\}\{{\boldsymbol{M}}\}$, and ${\rm vec}({\boldsymbol{M}})$, respectively. $\prod$ denotes the conventional or Cartesian product, depending on the factors. $\otimes$ denotes the Kronecker product. $\mathbbm{1}\{A\}$ denotes the indicator function whose value is 1 if $A$ is true and 0 otherwise. $[n]:=\{1,2,\dots,n\}$. $\propto$ means “proportional to”. The Grassmann manifold $G(\mathbb{C}^{T},K)$ is the space of $K$-dimensional subspaces in $\mathbb{C}^{T}$. In particular, $G(\mathbb{C}^{T},1)$ is the Grassmannian of lines. The Kullback-Leibler divergence of a distribution $p$ from another distribution $q$ of a random vector ${\boldsymbol{\mathsf{x}}}$ with domain ${\mathcal{X}}$ is defined by $D(q\|p):=\int_{{\mathcal{X}}}q({\boldsymbol{x}})\log\frac{q({\boldsymbol{x}})}{p({\boldsymbol{x}})}\mathop{}\mathrm{d}{\boldsymbol{x}}$ if ${\mathcal{X}}$ is continuous and $D(q\|p):=\sum_{{\boldsymbol{x}}\in{\mathcal{X}}}q({\boldsymbol{x}})\log\frac{q({\boldsymbol{x}})}{p({\boldsymbol{x}})}$ if ${\mathcal{X}}$ is discrete. ${\mathcal{N}}({\boldsymbol{\mu}},\boldsymbol{\Sigma})$ denotes the complex Gaussian vector distribution with mean ${\boldsymbol{\mu}}$, covariance matrix $\boldsymbol{\Sigma}$, and thus probability density function (PDF)

[TABLE]

II System Model

II-A Channel Model

We consider a SIMO MAC in which $K$ single-antenna users transmit to an $N$-antenna receiver. We assume that the channel is flat and block fading with an equal-length and synchronous (across the users) coherence interval of $T$ channel uses. That is, the channel vectors ${\boldsymbol{\mathsf{h}}}_{k}\in\mathbb{C}^{N\times 1}$, which contain the fading coefficients between the transmit antenna of user $k\in[K]$ and the $N$ receive antennas, remain constant within each coherence block of $T$ channel uses and change independently between blocks. Furthermore, the distribution of ${\boldsymbol{\mathsf{h}}}_{k}$ is assumed to be known to the receiver, but its realizations are unknown to both ends of the channel. Since the users are not co-located, we assume that the ${\boldsymbol{\mathsf{h}}}_{k}$ are independent across users. We consider Rayleigh fading with receiver-side correlation, i.e., ${\boldsymbol{\mathsf{h}}}_{k}\sim{\mathcal{N}}(\mathbf{0},\boldsymbol{\Xi}_{k})$, where $\boldsymbol{\Xi}_{k}\in\mathbb{C}^{N\times N}$ is the spatial correlation matrix. We assume that $\frac{1}{N}{\rm tr}\left\{\boldsymbol{\Xi}_{k}\right\}=:\xi_{k}$ where $\xi_{k}$ is the large-scale average channel gain from one of the receive antennas to user $k$. We assume that $T>K$ and $N\geq K$.

Within a coherence block, each transmitter $k$ sends a signal vector ${\boldsymbol{\mathsf{s}}}_{k}\in\mathbb{C}^{T}$, and the receiver receives a realization ${\boldsymbol{Y}}$ of the random matrix

[TABLE]

where ${\boldsymbol{\mathsf{S}}}=[{\boldsymbol{\mathsf{s}}}_{1}\ \dots\ {\boldsymbol{\mathsf{s}}}_{K}]\in\mathbb{C}^{T\times K}$ and ${\boldsymbol{\mathsf{H}}}=[{\boldsymbol{\mathsf{h}}}_{1}\ \dots\ {\boldsymbol{\mathsf{h}}}_{K}]\in\mathbb{C}^{N\times K}$ concatenate the transmitted signals and channel vectors, respectively, ${\boldsymbol{\mathsf{W}}}\in\mathbb{C}^{T\times N}$ is the Gaussian noise with i.i.d. ${\mathcal{N}}(0,\sigma^{2})$ entries independent of ${\boldsymbol{\mathsf{H}}}$, and the block index is omitted for simplicity.

We assume that the transmitted signals have average unit norm, i.e., $\mathbb{E}\left[\|{\boldsymbol{\mathsf{s}}}_{k}\|^{2}\right]=1,k\in[K]$. Under this normalization, the signal-to-noise ratio (SNR) of the transmitted signal from user $k$ at each receive antenna is ${\rm SNR}_{k}={\xi_{k}}/(T\sigma^{2})$. We assume that the transmitted signals belong to disjoint finite discrete individual constellations with vector-valued symbols. That is, ${\boldsymbol{\mathsf{s}}}_{k}\in{\mathcal{S}}_{k}:=\{{\boldsymbol{s}}_{k}^{(1)},\dots,{\boldsymbol{s}}_{k}^{(|{\mathcal{S}}_{k}|)}\},~{}k\in[K].$ In particular, ${\mathcal{S}}_{k}$ can be a Grassmannian constellation on $G(\mathbb{C}^{T},1)$, i.e., each constellation symbol ${\boldsymbol{s}}_{k}^{(i)}$ is a unit-norm vector representative of a point in $G(\mathbb{C}^{T},1)$. Another example is when the constellation symbols contain pilots and scalar data symbols.111In this case, the constellations are disjoint thanks to the fact that pilot sequences are user-specific. Each symbol in ${\mathcal{S}}_{k}$ is labeled with a binary sequence of length $B_{k}:=\log_{2}|{\mathcal{S}}_{k}|$.

II-B Multi-User Detection Problem

Given ${\boldsymbol{\mathsf{S}}}={\boldsymbol{S}}=[{\boldsymbol{s}}_{1},\ {\boldsymbol{s}}_{2},\ \dots,\ {\boldsymbol{s}}_{K}]$, the conditional probability density $p_{{\boldsymbol{\mathsf{Y}}}|{\boldsymbol{\mathsf{S}}}}$, also known as likelihood function, is derived similar to [27, Eq.(9)] as

[TABLE]

Given the received signal ${\boldsymbol{\mathsf{Y}}}={\boldsymbol{Y}}$, the joint multi-user ML symbol decoder is then

[TABLE]

Since the ML decoding metric depends on ${\boldsymbol{S}}$ only through $\sum_{k=1}^{K}{\boldsymbol{s}}_{k}{\boldsymbol{s}}_{k}^{\scriptscriptstyle\mathsf{H}}\otimes\boldsymbol{\Xi}_{k}$, for identifiability, it must hold that $\sum_{k=1}^{K}{\boldsymbol{s}}_{k}{\boldsymbol{s}}_{k}^{\scriptscriptstyle\mathsf{H}}\otimes\boldsymbol{\Xi}_{k}\neq\sum_{k=1}^{K}{\boldsymbol{s}}^{\prime}_{k}{{\boldsymbol{s}}^{\prime}_{k}}^{\scriptscriptstyle\mathsf{H}}\otimes\boldsymbol{\Xi}_{k}$ for any pair of distinct joint symbols ${\boldsymbol{S}}=[{\boldsymbol{s}}_{1},\dots,{\boldsymbol{s}}_{K}]$ and ${\boldsymbol{S}}^{\prime}=[{\boldsymbol{s}}^{\prime}_{1},\dots,{\boldsymbol{s}}^{\prime}_{K}]$ in $\prod_{k=1}^{K}{\mathcal{S}}_{k}$.

When a channel code is used, most channel decoders require the LLRs of the bits. The LLR of the $j$-th bit of user $k$, denoted by ${\mathsf{b}}_{k,j}$, given the observation ${\boldsymbol{\mathsf{Y}}}={\boldsymbol{Y}}$ is defined as

[TABLE]

where ${\mathcal{S}}_{k,j}^{(b)}$ denotes the set of all possible symbols in ${\mathcal{S}}_{k}$ with the $j$-th bit being equal to $b$ for $j\in[B_{k}]$ and $b\in\{0,1\}$. To compute (6), the posteriors $p_{{\boldsymbol{\mathsf{s}}}_{k}|{\boldsymbol{\mathsf{Y}}}}$, $k\in[K]$, are marginalized from

[TABLE]

Assuming that the transmitted signals are independent and uniformly distributed over the respective constellations, the prior $p_{{\boldsymbol{\mathsf{S}}}}$ factorizes as

[TABLE]

On the other hand, the likelihood function $p_{{\boldsymbol{\mathsf{Y}}}|{\boldsymbol{\mathsf{S}}}}({\boldsymbol{Y}}|[{\boldsymbol{s}}_{1},\dots,{\boldsymbol{s}}_{K}])$ involves all ${\boldsymbol{s}}_{1},\dots,{\boldsymbol{s}}_{K}$ in such a manner that it does not straightforwardly factorize. Exact marginalization of $p_{{\boldsymbol{\mathsf{S}}}|{\boldsymbol{\mathsf{Y}}}}$ requires computing

[TABLE]

That is, it requires computing $p_{{\boldsymbol{\mathsf{Y}}}|{\boldsymbol{\mathsf{S}}}}({\boldsymbol{Y}}|{\boldsymbol{S}})$ (which requires the inversion of an $NT\times NT$ matrix) for all ${\boldsymbol{S}}\in\prod_{k=1}^{K}{\mathcal{S}}_{k}$. Thus, the total complexity of exact marginalization is $O({K^{6}}2^{KB})$.222Throughout the paper, as far as the complexity analysis is concerned, we assume for notational simplicity that $T=O(K)$, $N=O(K)$, and $|{\mathcal{S}}_{k}|=O(2^{B}),~{}\forall k\in[K]$. If the channels are uncorrelated ($\boldsymbol{\Xi}_{k}={\boldsymbol{I}}_{N}$), the likelihood function can be simplified as $p_{{\boldsymbol{\mathsf{Y}}}|{\boldsymbol{\mathsf{S}}}}({\boldsymbol{Y}}|{\boldsymbol{S}})=\frac{\exp\left(-{\rm tr}\left\{\right\}\left\{{\boldsymbol{Y}}^{\scriptscriptstyle\mathsf{H}}(\sigma^{2}{\boldsymbol{I}}_{T}+{\boldsymbol{S}}{\boldsymbol{S}}^{\scriptscriptstyle\mathsf{H}})^{-1}{\boldsymbol{Y}}\right\}\right)}{\pi^{NT}{\rm det}\left(\right)^{N}(\sigma^{2}{\boldsymbol{I}}_{T}+{\boldsymbol{S}}{\boldsymbol{S}}^{\scriptscriptstyle\mathsf{H}})}$. Thus, the complexity of exact marginalization is reduced to $O(K^{3}2^{KB})$. This is formidable for many users or large constellations. Thus, we seek alternative approaches to estimate

[TABLE]

II-C Achievable Rate

According to [28, Sec. II], the highest sum-rate reliably achievable with a given decoding metric $\hat{p}_{{\boldsymbol{\mathsf{S}}}|{\boldsymbol{\mathsf{Y}}}}$, so-called the mismatched sum-rate, is lower bounded by the generalized mutual information (GMI) given by *

$t_{\max}$ denotes the number of iterations. $t_{0}\in[t_{\max}]$. *

VI-B Stabilization

We discuss some possible numerical problems in the EP algorithm and our solutions.

VI-B1 Singularity of $\boldsymbol{\Sigma}_{k}$

First, in (LABEL:eq:Sigma_k), since the $NT\times NT$ matrix $\boldsymbol{\Sigma}_{k}$ is the weighted sum of the terms of rank less than $NT$, it can be close to singular if at a certain iteration, only few of the weights $\pi_{k1}^{(i)}$ are sufficiently larger than zero. The singularity of $\boldsymbol{\Sigma}_{k}$ can also arise from the constellation structure. For example, the constellations proposed in **[22*]** are precoded versions of a constellation in $G(\mathbb{C}^{T-K+1},1)$ and the maximal rank of $\boldsymbol{\Sigma}_{k}$ is $N(T-K+1)\leq NT$. To avoid the inverse of $\boldsymbol{\Sigma}_{k}$, we express ${\boldsymbol{C}}_{k1}$ in (LABEL:eq:C_k1) and ${\boldsymbol{\mu}}_{k1}$ in (LABEL:eq:mu_k1) respectively as *

[TABLE]

VI-B2 “Negative variance”

Another problem is that ${\boldsymbol{C}}_{k1}$ is not guaranteed to be positive definite even if both ${\boldsymbol{C}}_{k0}$ and $\boldsymbol{\Sigma}_{k}$ are. When ${\boldsymbol{C}}_{k1}$ is not positive definite, from (LABEL:eq:C_k0), ${\boldsymbol{C}}_{k0}$ can have negative eigenvalues, which, through (LABEL:eq:pihat_k1), can make $\hat{\pi}_{k1}^{({i}_{k})}$ become close to a Kronecker-delta distribution (even at low SNR) where the position of the mode can be arbitrary, and the algorithm may diverge. Note that this “negative variance” problem is common in EP (see, e.g., **[24, Sec. 3.2.1]**, **[31, Sec. 5.3]**). There has been no generally accepted solution and one normally resorts to various heuristics adapted to each problem. In our problem, to control the eigenvalues of ${\boldsymbol{C}}_{k1}$, we modify (94) by first computing the eigendecomposition -{\boldsymbol{C}}_{k0}\big{(}\boldsymbol{\Sigma}_{k}-{\boldsymbol{C}}_{k0}\big{)}^{-1}\boldsymbol{\Sigma}_{k}={\boldsymbol{V}}\boldsymbol{\Lambda}{\boldsymbol{V}}^{-1}, then computing ${\boldsymbol{C}}_{k1}$ as ${\boldsymbol{C}}_{k1}={\boldsymbol{V}}|\boldsymbol{\Lambda}|{\boldsymbol{V}}^{-1},$ where $|\boldsymbol{\Lambda}|$ is the element-wise absolute value of $\boldsymbol{\Lambda}$. This manipulation of replacing the variance parameters by their absolute values was also used in **[32]**.

VI-B3 Overconfidence at early iterations

Finally, due to the nature of the message passing between continuous and discrete distribution, it can happen that all the mass of the PMF $\hat{\pi}_{k1}^{({i}_{k})}$ is concentrated on a small region of a potentially large constellation ${\mathcal{S}}_{k}$. For example, if ${\pi}_{k1}^{({i}_{k})}$ is close to a Kronecker-delta distribution with a single mode at $i_{0}$, then (LABEL:eq:xhat_ki) and (LABEL:eq:Sigma_ki) implies that $\boldsymbol{\Sigma}_{k}$ is approximately $\boldsymbol{\Sigma}_{ki_{0}}$, and then from (LABEL:eq:C_k1), ${\boldsymbol{C}}_{k1}\approx({\boldsymbol{s}}_{k}^{(i_{0})}{\boldsymbol{s}}_{k}^{(i_{0}){\scriptscriptstyle\mathsf{H}}})\otimes{\boldsymbol{\Xi}_{k}}$. In this case, almost absolute certainty is placed on the symbol ${\boldsymbol{s}}_{k}^{(i_{0})}$, and the algorithm will not be able significantly update its belief in the subsequent iterations. This can be problematic when the mode of ${\pi}_{k1}^{({i}_{k})}$ is placed on the wrong symbol at early iterations. To smooth the updates, we apply damping on the update of the parameters of the continuous distributions ${\mathcal{N}}({\boldsymbol{z}}_{k};{\boldsymbol{\mu}}_{k1},{\boldsymbol{C}}_{k1})$ and ${\mathcal{N}}({\boldsymbol{z}}_{k};{\boldsymbol{\mu}}_{k0},{\boldsymbol{C}}_{k0})$. That is, with a damping factor $\eta\in[0;1]$, at iteration $t$ and for each user $k$, we update

[TABLE]

*In short, we stabilize the EP message updates by replacing (96), (97), (98), and (99) for (LABEL:eq:C_k1), (LABEL:eq:mu_k1), (LABEL:eq:C_k0), and (LABEL:eq:mu_k0), respectively. This technique also applies to EPAK. For MMSE-SIA, we damp the update of ${\boldsymbol{Q}}_{k}$ and ${\boldsymbol{R}}_{k}$ in a similar manner as {\boldsymbol{Q}}_{k}(t)=\eta\big{(}\sum_{l\neq k}{\boldsymbol{R}}_{l}(t-1)+\sigma^{2}{\boldsymbol{I}}_{T}\big{)}+(1-\eta){\boldsymbol{Q}}_{k}(t-1) and ${\boldsymbol{R}}_{k}(t)=\eta\sum_{i_{k}=1}^{|{\mathcal{S}}_{k}|}\pi_{k1}^{(i_{k})}(t){\boldsymbol{s}}_{k}^{(i_{k})}{\boldsymbol{s}}_{k}^{(i_{k}){\scriptscriptstyle\mathsf{H}}}+(1-\eta){\boldsymbol{R}}_{k}(t-1).$ Note that damping does not change the complexity order of these schemes. The approaches described in this subsection were implemented for the numerical results in the next section. *

VII Performance Evaluation

In this section, we evaluate the performance of our proposed schemes for a given set of individual constellations. We assume that $B_{1}=\dots B_{K}=:B$. We consider the local scattering model **[4, Sec. 2.6]** for the correlation matrices $\boldsymbol{\Xi}_{k}$. Specifically, the $(l,m)$-th element of $\boldsymbol{\Xi}_{k}$ is generated as $[\boldsymbol{\Xi}_{k}]_{l,m}=\xi_{k}\mathbb{E}\left[\right]_{\delta_{k}}[\exp(2\pi d_{H}(l-m)\sin(\varphi_{k}+\delta_{k}))],$ where $d_{H}$ is the antenna spacing in the receiver array (measured in number of wavelengths), $\varphi_{k}$ is a deterministic nominal angle, and $\delta_{k}$ is a random deviation. We consider $d_{H}=\frac{1}{2}$, $\varphi_{k}$ generated uniformly in $[-\pi,\pi]$, and $\delta_{k}$ uniformly distributed in $[-\sqrt{3}\sigma_{\varphi},\sqrt{3}\sigma_{\varphi}]$ with angular standard deviation $\sigma_{\varphi}=10^{\circ}$. We also consider $\xi_{k}=1,\forall k.$ We set a damping factor $\eta=0.9$ for EP, EPAK, and MMSE-SIA.

VII-A Test Constellations, State-of-the-Art Detectors, and Benchmarks

VII-A1 Precoding-Based Grassmannian Constellations

We consider the constellation design in **[22]**, which imposes a geometric separation between the individual constellations through a set of precoders ${\boldsymbol{U}}_{k}$, $k\in[K]$. Specifically, starting with a Grassmannian constellation {\mathcal{D}}=\big{\{}{\boldsymbol{d}}^{(1)},\dots,{\boldsymbol{d}}^{(2^{B})}\} in $G(\mathbb{C}^{T-K+1},1)$, the individual constellation ${\mathcal{S}}_{k}$ is generated as

[TABLE]

We consider the precoders ${\boldsymbol{U}}_{k}$ defined in **[22, Eq.(11)]** and two candidates for ${\mathcal{D}}$:

• •

A numerically optimized constellation generated by solving the max-min distance criteria

[TABLE]

where $d({\boldsymbol{d}}^{(i)},{\boldsymbol{d}}^{(j)}):=\sqrt{1-|{\boldsymbol{d}}^{(i){\scriptscriptstyle\mathsf{H}}}{\boldsymbol{d}}^{(j)}|^{2}}$ is the chordal distance between two Grassmannian points represented by ${\boldsymbol{d}}^{(i)}$ and ${\boldsymbol{d}}^{(j)}$. A constellation with maximal minimum pairwise distance leads to low symbol error rate in the absence of the interference. In our simulation, we approximate (101) by \min_{{\mathcal{D}}}\log\sum_{1\leq i<j\leq 2^{B}}\exp\big{(}\frac{|{\boldsymbol{d}}^{(i){\scriptscriptstyle\mathsf{H}}}{\boldsymbol{d}}^{(j)}|}{\epsilon}\big{)} with a small $\epsilon$ for smoothness, then solve it using gradient descent on the Grassmann manifold using the Manopt toolbox **[33]**.

• •

The cube-split constellation proposed in **[34, 17]**. This structured constellation has good distance properties and allows for low-complexity single-user decoding and a simple yet effective binary labeling scheme.

Exploiting the precoder structure, **[22]** introduced a detector **[22, Sec. V-B-3]** that iteratively mitigates interference by projecting the received signal onto the subspace orthogonal to the interference subspace. We refer to it as POCIS (Projection onto the Orthogonal Complement of the Interference Subspace). For each user $k$, POCIS first estimates the row space of the interference $\sum_{l\neq k}{\boldsymbol{\mathsf{s}}}_{l}{\boldsymbol{\mathsf{h}}}_{l}^{\scriptscriptstyle\mathsf{T}}$ based on the precoders and projects the received signal onto the orthogonal complement of this space. It then performs single-user detections to obtain point estimates of the transmitted symbols. From these estimates, POCIS estimates the column space of the interference and projects the received signal onto its orthogonal complement. This process is repeated in the next iteration. The complexity order of POCIS is equivalent to the MMSE-SIA scheme. Note that only the indices of the estimated symbols are passed in POCIS, as opposed to the soft information on the symbols as in EP, MMSE-SIA, and EPAK.

VII-A2 Pilot-Based Constellations

We also consider the pilot-based constellations in which the symbols are generated as {\boldsymbol{s}}_{k}^{(i)}=\Big{[}\sqrt{\frac{K}{T}}{\boldsymbol{e}}_{k}^{\scriptscriptstyle\mathsf{T}}\ \sqrt{\frac{T-K}{TP_{\rm avg}}}\tilde{{\boldsymbol{s}}}_{k}^{(i){\scriptscriptstyle\mathsf{T}}}\Big{]}^{\scriptscriptstyle\mathsf{T}} where ${\boldsymbol{e}}_{k}$ is the $k$-th column of ${\boldsymbol{I}}_{K}$, $\tilde{{\boldsymbol{s}}}_{k}^{(i)}$ is a vector of data symbols taken from a scalar constellation, such as QAM, and $P_{\rm avg}$ is the average symbol power of the considered scalar constellation. Note that this corresponds to the scenario where the $K$ users transmit mutually orthogonal pilot sequences, followed by spatially multiplexed parallel data transmission. Many MIMO detectors have been proposed specifically for these constellations. We consider some representatives as follows.

• •

The receiver MMSE-estimates the channel based on the first $K$ rows of ${\boldsymbol{\mathsf{Y}}}$, then MMSE-equalizes the received data symbols in the remaining $T-K$ rows of ${\boldsymbol{\mathsf{Y}}}$, and performs a scalar demapper on the equalized symbols.

• •

The receiver MMSE-estimates the channel, then decodes the data symbols using the Schnorr-Euchner sphere decoder **[35]**, referred to as SESD.

• •

The receiver performs the semi-blind joint ML channel estimation and data detection scheme in **[9]** with repeated weighted boosting search (RWBS) for channel estimation and the Schnorr-Euchner sphere decoder for data detection, referred to as RWBS-SESD.

We note that the sphere decoder has near optimal performance given the channel knowledge, but its complexity is non-deterministic and can be exponential in the channel dimension if the channel matrix is ill-conditioned.

VII-A3 Benchmarks

We consider the optimal ML detector, whenever it is feasible, as a benchmark. When the optimal detector is computationally infeasible, we resort to another benchmark consisting in giving the receiver, while it decodes the signal ${\boldsymbol{\mathsf{s}}}_{k}$ of user $k$, the knowledge of the signals ${\boldsymbol{\mathsf{s}}}_{l}$ (but not the channel ${\boldsymbol{\mathsf{h}}}_{l}$) of all the interfering users $l\neq k$. With this genie-aided information, optimal ML decoding (LABEL:eq:MLdecoder) can be performed by keeping ${\boldsymbol{\mathsf{s}}}_{l}$ fixed for all $l\neq k$ and searching for the best ${\boldsymbol{\mathsf{s}}}_{k}$ in ${\mathcal{S}}_{k}$, thus reducing the total search space size from $2^{BK}$ to $K2^{B}$. The posterior marginals are computed separately for each user accordingly. This genie-aided detector gives an upper bound on the performance of EP, MMSE-SIA, EPAK, and POCIS.

VII-B Convergence and Running Time

To assess the convergence of the algorithms, we evaluate the total variation distance between the estimated marginal posteriors $\hat{p}_{{\boldsymbol{\mathsf{s}}}_{k}|{\boldsymbol{\mathsf{Y}}}}$ at each iteration and the exact marginal posteriors $p_{{\boldsymbol{\mathsf{s}}}_{k}|{\boldsymbol{\mathsf{Y}}}}$ when exact marginalization (9) is possible. The total variation distance between two probability measures $P$ and $Q$ on ${\mathcal{X}}$ is defined as ${\mathcal{T}}{\mathcal{V}}(P,Q):=\frac{1}{2}\sum_{x\in{\mathcal{X}}}|P(x)-Q(x)|$. At iteration $t$ where the estimated posteriors are $\hat{p}^{(t)}_{{\boldsymbol{\mathsf{s}}}_{k}|{\boldsymbol{\mathsf{Y}}}}$, $k\in[K]$, we evaluate the average total variation distance as

[TABLE]

We consider the precoding-based Grassmannian constellations. Fig. 3 shows the empirical average total variation $\Delta_{t}$ for $T=6$, $K=3$, $N=4$, and $B=4$ at ${\rm SNR}=8$ dB. As can be seen, at convergence, EP provides the most accurate estimates of the marginal posteriors although it is less stable than other schemes. EP converges after $6$ iterations while MMSE-SIA converges after $5$ iterations. For uncorrelated fading, EPAK with $t_{0}=2$ can be eventually better than MMSE-SIA, but converges slower. For correlated fading, EPAK totally fails because of the inaccuracy of the approximation with Kronecker products. POCIS converges very quickly after $2$ iterations but achieves a relatively low accuracy of the posterior estimation.

Fig. 4 depicts the average running time (on a local server) of exact marginalization compared with $6$ iterations of EP, EPAK, MMSE-SIA, and POCIS at ${\rm SNR}=8$ dB. These schemes have significantly lower computation time than exact marginalization. The running time saving of EPAK w.r.t. EP is not significant, even with $t_{0}=0$. For uncorrelated fading, MMSE-SIA has much shorter running time than all other schemes.

From these convergence behaviors, hereafter, we fix the number of iterations of EP, MMSE-SIA, and EPAK as $6$ and of POCIS as $3$. Furthermore, we consider EPAK only for uncorrelated fading. For correlated fading, we generate the correlation matrices once and fix them over the simulation.

VII-C Achievable Rate

We first plot the achievable mismatched sum-rate $R_{\rm GMI}$ of the system calculated as in (LABEL:eq:approx_achievableRate) for $T=6$, $K=3$, $N=4$ and $B\in\{4,8\}$ in Fig. 5. We consider the precoding-based Grassmannian constellations. For ${\mathcal{D}}$, we use the numerically optimized constellation if $B=4$ and the cube-split constellation if $B=8$. For uncorrelated fading (Fig. 5(a), the rates achieved with EP and MMSE-SIA detectors are very close to the achievable rate of the system (with the optimal detector) and not far from that of the genie-aided detector. EPAK (with $t_{0}=2$) achieves a very low rate, especially in the low SNR regime where the Kronecker approximation is not accurate. For correlated fading, (Fig. 5(b)), the rates achieved with EP and MMSE-SIA are only marginally lower than that of the optimal detector and genie-aided detector. In both cases, the rate achieved with POCIS is lower than that of EP and MMSE-SIA in the lower SNR regime and converges slowly with SNR to the limit $\frac{BK}{T}$ bits/channel use.

VII-D Symbol Error Rates of Hard Detection

*Next, we use the outputs of EP, EPAK, MMSE-SIA and POCIS for a maximum-a-posteriori (MAP) hard detection. We evaluate the performance in terms of symbol error rate (SER). *

In Fig. 6, we consider the precoding-based constellations with $T=6$, $K=3$, $N\in\{4,8\}$, and $B=4$, for which the optimal ML detector (LABEL:eq:MLdecoder) is computationally feasible. We observe that the SER of the EP and MMSE-SIA detectors are not much higher than that of the optimal detector, especially in the lower SNR regime. The SER of EPAK is significantly higher than that of EP and MMSE-SIA for $t_{0}=0$. This is greatly improved by setting $t_{0}=2$, i.e., keeping the first two iterations of EP. The gain of EP w.r.t. EPAK and MMSE-SIA is more pronounced when the SNR increases. For correlated fading, EP performs almost as good as the optimal detector, whose SER performance is closely approximated by the genie-aided detector.

In Fig. 7, we consider $T=6$, $K=3$, $N=8$, and $B=9$ and use the genie-aided detector as a benchmark. In Fig. 7(a), we consider uncorrelated fading and use the pilot-based constellations with $8$-QAM data symbols. The performance of EP is very close to that of the genie-aided detector. The performance of MMSE-SIA is close to EP in the low SNR regime (${\rm SNR}\leq 8$* dB). We also depict the SER of the three pilot-based detectors in Section VII-A2, namely, 1) MMSE channel estimation, MMSE equalizer, and QAM demapper, 2) SESD, and 3) RWBS-SESD. These three schemes are outperformed by the EP detectors. In Fig. 7(b), we consider correlated fading and use the precoding-based Grassmannian constellations with ${\mathcal{D}}$ numerically optimized. We observe again that EP achieves almost the same SER performance as the genie-aided detector.*

VII-E Bit Error Rates with a Channel Code

In this subsection, we use the output of the soft detectors for channel decoding. We consider the precoding-based Grassmannian constellations with the cube-split constellation for ${\mathcal{D}}$ since it admits an effective and simple binary labeling **[17]**. We take the binary labels of the symbols in ${\mathcal{D}}$ for the corresponding symbols in ${\mathcal{S}}_{k}$. We integrate a standard symmetric parallel concatenated rate-$1/3$* turbo code [36]. The turbo encoder accepts packets of $1008$ bits; the turbo decoder computes the bit-wise LLR from the soft outputs of the detection scheme as in (6) and performs 10 decoding iterations for each packet.*

In Fig. 8, we show the bit error rate (BER) with this turbo code using $B=8$ bits/symbol and different values of $T$ and $K=N$. EP achieves the closest performance to the genie-aided detector and the optimal detector (9). The BER of MMSE-SIA vanishes slower with the SNR than the other schemes, and becomes better than POCIS as $K$ and $N$ increase. The BER of EPAK with $t_{0}=2$ is higher than all other schemes. Under uncorrelated fading, for $T=7$ and $K=N=4$, the power gain of EP w.r.t. MMSE-SIA, POCIS, and EPAK for the same BER of $10^{-3}$ is about $3$ dB, $4$ dB, and $8$ dB, respectively. We also observe that the genie-aided detector gives very optimistic BER performance results compared to the optimal detector.

Finally, in Fig. 9, we consider $T=6$, $K=3$, $N=4$, and compare the BER with the same turbo code for different $B$. For $B=5$, both EP and MMSE-SIA have performance close to the optimal detector. Under uncorrelated fading, MMSE-SIA can be slightly better than EP. This is due to the residual effect (after damping) of the phenomenon that all the mass of $\pi_{k1}^{({i}_{k})}$ is concentrated on a possibly wrong symbol at early iterations, and EP may not be able to refine significantly the PMF in the subsequent iterations if the constellation is sparse. This situation is not observed for $B=8$, i.e., larger constellations. Also, as compared to the case $T=6,K=3,B=8$ in Fig. 8, the performance of MMSE-SIA is significantly improved as the number of receive antennas increases from $N=3$ to $N=4$. As in the previous case, EPAK does not perform well.

VIII Conclusion

We proposed an expectation propagation (EP) based scheme and two simplifications (EPAK and MMSE-SIA) of this scheme for multi-user detection in non-coherent SIMO multiple access channel with spatially correlated Rayleigh fading. EP and MMSE-SIA are shown to achieve good performance in terms of mismatched sum-rate, symbol error rate when they are used for hard detection, and bit error rate when they are used for soft-input soft-output channel decoding. EPAK has acceptable performance with uncorrelated fading. It performs well for hard symbol detection but inadequately for soft-output detection. While MMSE-SIA and EPAK have lower complexity than EP, the performance gain of EP with respect to MMSE-SIA and EPAK is more significant when the number of users and/or the constellation size increase. Possible extensions of this work include considering more complicated fading models and analyzing theoretically the performance of EP for non-coherent reception.

-A Properties of the Gaussian PDF

Lemma 1**.**

Let ${\boldsymbol{\mathsf{x}}}$ be an $n$-dimensional complex Gaussian vector. It holds that

${\mathcal{N}}({\boldsymbol{x}};{\boldsymbol{\mu}},\boldsymbol{\Sigma})={\mathcal{N}}({\boldsymbol{x}}+{\boldsymbol{y}};{\boldsymbol{\mu}}-{\boldsymbol{y}},\boldsymbol{\Sigma})$* for ${\boldsymbol{y}}\in\mathbb{C}^{n}$;* 2. 2.

Gaussian PDF multiplication rule: ${\mathcal{N}}({\boldsymbol{x}};{\boldsymbol{\mu}}_{1},\boldsymbol{\Sigma}_{1}){\mathcal{N}}({\boldsymbol{x}};{\boldsymbol{\mu}}_{2},\boldsymbol{\Sigma}_{2})={\mathcal{N}}({\boldsymbol{x}};{\boldsymbol{\mu}}_{{\rm new}},\boldsymbol{\Sigma}_{{\rm new}}){\mathcal{N}}(\mathbf{0};{\boldsymbol{\mu}}_{1}-{\boldsymbol{\mu}}_{2},\boldsymbol{\Sigma}_{1}+\boldsymbol{\Sigma}_{2})$, where \boldsymbol{\Sigma}_{{\rm new}}:=\big{(}\boldsymbol{\Sigma}_{1}^{-1}+\boldsymbol{\Sigma}_{2}^{-1}\big{)}^{-1} and {\boldsymbol{\mu}}_{{\rm new}}:=\boldsymbol{\Sigma}_{{\rm new}}\big{(}\boldsymbol{\Sigma}_{1}^{-1}{\boldsymbol{\mu}}_{1}+\boldsymbol{\Sigma}_{2}^{-1}{\boldsymbol{\mu}}_{2}\big{)}.

Proof.

The first part follows readily from the definition of ${\mathcal{N}}({\boldsymbol{x}};{\boldsymbol{\mu}},\boldsymbol{\Sigma})$. The complex Gaussian PDF multiplication rule is a straightforward generalization of the real counterpart [37].

∎

-B Proof of Proposition LABEL:prop:pnew

Using the natural logarithm for the KL divergence, we derive

[TABLE]

[TABLE]

where (111) follows from q_{\alpha}({\boldsymbol{x}})=\frac{\psi_{\alpha}({\boldsymbol{x}}_{\alpha})}{m_{\alpha}({\boldsymbol{x}}_{\alpha})}\Big{[}\prod_{\beta\in\mathfrak{N}_{\alpha}}\hat{p}_{\beta}({\boldsymbol{x}}_{\beta})\Big{]}\Big{[}\prod_{\beta\notin\mathfrak{N}_{\alpha}}\hat{p}_{\beta}({\boldsymbol{x}}_{\beta})\Big{]}, and (112) follows from (LABEL:eq:expfam). From (LABEL:eq:KLalf1), we can see that the optimization (LABEL:eq:pnew) of $\underline{p}$ decouples over $\underline{p}_{\beta}$, and the optimal distribution can be expressed as $\hat{p}^{{\rm new}}_{\alpha}({\boldsymbol{x}})=\prod_{\beta}\hat{p}^{{\rm new}}_{\alpha,\beta}({\boldsymbol{x}}_{\beta})$. For $\beta\notin\mathfrak{N}_{\alpha}$, the minimum of D\big{(}\hat{p}_{\beta}\big{\|}\underline{p}_{\beta}\big{)} is simply [math] and achieved with $\hat{p}^{{\rm new}}_{\alpha,\beta}({\boldsymbol{x}}_{\beta})={\hat{p}}_{\beta}({\boldsymbol{x}}_{\beta})$. For $\beta\in\mathfrak{N}_{\alpha}$, since the log-partition function $A_{\beta}(\underline{{\boldsymbol{\gamma}}}_{\beta})$ is convex in $\underline{{\boldsymbol{\gamma}}}_{\beta}$ (see, e.g., **[38, Lemma 1]**), the minimum of A_{\beta}(\underline{{\boldsymbol{\gamma}}}_{\beta})-\underline{{\boldsymbol{\gamma}}}_{\beta}^{\scriptscriptstyle\mathsf{T}}\mathbb{E}\left[\right]_{q_{\alpha}}\big{[}{\boldsymbol{\phi}}({\boldsymbol{x}}_{\beta})\big{]} is achieved at the value of $\underline{\bm{{\boldsymbol{\gamma}}}}_{\beta}$ where its gradient is zero. Using the well-known property of the log-partition function, $\nabla_{{\boldsymbol{\gamma}}}A_{\beta}({\boldsymbol{\gamma}})=\mathbb{E}\left[\right]_{\hat{p}_{\beta}}[{\boldsymbol{\phi}}_{\beta}({\boldsymbol{\gamma}})]$, we get that the zero-gradient equation is equivalent to the moment matching criterion $\mathbb{E}\left[\right]_{\hat{p}^{{\rm new}}_{\alpha,\beta}}[{\boldsymbol{\phi}}_{\beta}({\boldsymbol{x}}_{\beta})=\mathbb{E}\left[\right]_{q_{\alpha}}[{\boldsymbol{\phi}}_{\beta}({\boldsymbol{x}}_{\beta})]$.