1 Introduction↩︎

Massive multiple input multiple output (mMIMO) [1] is a key technology for modern wireless communications. By equipping mMIMO base stations (BS) with a large number of antennas, it can handle a large number of users that transmit on the same time/frequency resource; the increased spatial multiplexing brings massive gain in spectral and energy efficiency. However, this entails joint detection of transmitted symbols from a high-dimensional and sometimes ill-conditioned channel. This high dimensionality requires algorithms capable of achieving a satisfactory trade-off between performance and complexity.

The optimal maximum a posteriori (MAP) detection rapidly becomes untractable because it has exponential complexity in constellation size and the number of symbols transmitted. Most of the algorithms introduced in the literature overcome this issue by assuming that received symbols follow a Gaussian distribution; this relaxation leads e.g. to zero forcing (ZF) or linear minimum mean square error (L-MMSE) detection, which however is still computationally intensive for large systems since it requires the computation of the pseudo-inverse of the channel matrix. Gaussian relaxation has also been used to define approximate iterative algorithms schemes such as Gaussian belief propagation (Ga-BP) [2], expectation propagation (EP) [3] or approximate message passing (AMP) [4]. A good introduction to the application of these message passing algorithms to mMIMO detection and the relation between these three algorithms can be found in [5].

In this paper, we propose a new iterative detection algorithm applicable to PSK-modulated mMIMO based on directional statistics [6], more specifically von Mises (vM) distributions. The use of directional statistics for the physical layer was studied in [7] to decode a coded pulse amplitude modulation. Building on variational message passing [8] and mean field BP [9], [10], a BP algorithm using von Mises distributions was introduced in [11] in order to decode tensor-based modulations for unsourced random access.

The paper is organized as follows. We first introduce the detection problem in mMIMO in section 2. In section 3, we briefly review the properties of von Mises distributions, and derive a message-passing algorithm that leverage von Mises belief propagation (vM-BP), for coherent detection of PSK symbols for the case of perfect channel state information (CSI). Section 4 addresses the case of imperfect CSI. In section 5, we detail the performance and complexity of the proposed approach, and benchmark it against state-of-the-art methods. Finally, section 6 concludes the article.

2 System Model and Problem Formulation↩︎

We consider a MIMO system with \(K\) single-antenna transmitters and a receiver equipped with \(N\) antennas. At each channel access, every user transmits \(J\) bits of information that are mapped to a discrete PSK constellation \(\Omega = \{e^{j2\pi m/M}, m=0,\dots,M-1 \}\) of cardinality \(M=2^J\) (M-PSK). For simplicity and without loss of generality, we assume that every user uses the same constellation. The symbol transmitted by user \(k\) is denoted by \(x^{(k)}, k = 1,\ldots,K\). We assume a narrow-band fading channel, with gains denoted by \(\mathbf{H} = [\mathbf{h}^{(1)},\dots ,\mathbf{h}^{(K)}]\in\mathbb{C}^{N\times K}\) where \(\mathbf{h}^{(k)} \in \mathbb{C}^N\) is the channel gains between user \(k\) and \(N\) receive antennas. The received signal is hence given by \[\mathbf{y} = \sum_{k=1}^{K} x^{(k)}\mathbf{h}^{(k)} + \mathbf{w}, \label{MIMO95rx95signal}\tag{1}\] where \(\mathbf{y} \in \mathbb{C}^{N}\) and \(\mathbf{w} \in \mathbb{C}^N\) is a complex circular additive white Gaussian noise (AWGN) with zero mean and covariance matrix \(\sigma^2 \mathbf{I}_N\), where \(\mathbf{I}_N\) denotes the identity matrix of size \(N\).

Assuming that \(\mathbf{H}\) is known, the optimal solution for mMIMO detection considering the MAP criterion is given by: \[\label{eq95MAP} \hat{\mathbf{x}} = \arg\max_{\mathbf{x} \in \Omega^{K}} p(\mathbf{x} | \mathbf{y}; \mathbf{H} ).\tag{2}\] We can write the corresponding vector factorization of a posteriori vector probability as \[\begin{align} \label{eq:vector95factor} p(\mathbf{x} | \mathbf{y}; \mathbf{H}) & \propto p(\mathbf{y} | \mathbf{x};\mathbf{H}) p(\mathbf{x}) \nonumber\\ & \propto p(\mathbf{y} | \mathbf{x};\mathbf{H}) \prod_{k=1}^{K} f_{X^{(k)}}(x^{(k)}), \end{align}\tag{3}\] where \(p(\mathbf{y} | \mathbf{x}; \mathbf{H})\) is the likelihood associated with the observation model given by 1 and \(f_{X^{(k)}}(x)\) is the prior distribution for variable \(X^{(k)}\) (usually assumed uniform over \(\Omega\)). Unfortunately, computing \(\hat{\mathbf{x}}\) as per 2 has prohibitive complexity in mMIMO systems, because of the cardinality of the search space given by \(\left|\Omega^{K}\right| = 2^{JK}\), which scales exponentially with \(K\) and \(J\). To overcome this complexity, suboptimal criteria are classically considered based on linear block detection, such as L-MMSE. Alternatively, approximate message-passing based detection methods are often considered to iteratively compute approximations of the a posteriori marginals \(p\left(x^{(k)} \mid \mathbf{y};\mathbf{H}\right)\), enabling approximate MAP symbol detection. In the MIMO context, capitalizing on the knowledge of the channel matrix \(\mathbf{H}\), we can rewrite 3 as \[\begin{align} p(\mathbf{x} \mid \mathbf{y}; \mathbf{H}) & \propto \prod_{n=1}^{N} p\left(y_n\mid \mathbf{x}; \mathbf{H}\right) \prod_{k=1}^{K} f_{X^{(k)}}(x^{(k)}) \nonumber \\ & \propto \prod_{n=1}^{N} f_n(\mathbf{x}; y_n,\mathbf{h}_n) \prod_{k=1}^{K} f_{X^{(k)}}(x^{(k)}), \label{eq:scalar95factor} \end{align}\tag{4}\] where \(f_n(\mathbf{x}; y_n,\mathbf{h}_n)\triangleq p\left(y_n \mid \mathbf{x}; \mathbf{H}\right) = p\left(y_n \mid \mathbf{x}; \mathbf{h}_n\right)\), with \(\mathbf{h}_n\) being the \(n\)-th row of \(\mathbf{H}\). This factorization naturally lends itself to a bipartite factor graph representation [12], for which latent variable nodes are associated with symbols \(x^{(k)}\) and observed variable nodes are associated with antenna observations \(y_n\). The receive antenna \(n\) is associated with a factor node \(f_{n}\) (also denoted function node), defined by the corresponding local likelihood. This factor node links the observed variable \(y_n\) to the latent variables \(x^{(k)}, k=1, \ldots, K\), that contribute to this observation through the channel.

Based on the factor graph representation, different message-passing strategies can be considered to further lower the computational complexity. A common approach is to constrain the exchanged messages, which consist in distributions describing the information about a variable aggregated from the neighboring nodes in the graph, to belong to the family of Gaussian distributions. Doing so, each message is entirely characterized by a small number of parameters, such as its mean and variance. Therefore, message-passing no longer consists of transmitting complete distributions, but only their parameters. Such parametric message-passing approaches on the factor graph associated to the factorization in 4 lead to scalar algorithms such as Ga-BP, EP, and AMP. Note that when we consider the vector factorization as in 3 , applying message-passing rule updates on the corresponding factor graph leads to the vector counterpart of the previously cited algorithm, such as the vector EP (VEP) algorithm [13].

In this paper, we restrict our study to M-PSK modulations, for which all symbols lie on the unit circle. This motivates the use of directional instead of Gaussian distributions for the parametric messages, as explained next.

3 vM-BP detection for Coherent MIMO↩︎

3.1 von Mises Distribution: A Brief Review↩︎

The von Mises distribution (also known as circular normal distribution) [6] is a distribution defined on the unit circle \(\mathcal{C} = \{x \in \mathbb{C} \;\mathrm{s.t.} \;|x|=1 \}\) which belongs to the exponential family. Let \(X \in \mathcal{C}\) be a random variable following a von Mises distribution with parameter \(\eta \in \mathbb{C}\); by definition, the associated probability density function (p.d.f.) is \[\mathrm{vM}(x;\eta)=\frac{1}{2 \pi I_0(\lvert \eta \rvert)} \exp\big(\mathrm{Re}(\eta^*x)\big).\] \(\eta\) can be written as \(\eta=\kappa e^{j\mu}\), where \(\mu = \mathrm{arg}(\eta) \in [0,2\pi)\) and \(\kappa = \lvert \eta \rvert \in [0,\infty)\) are respectively the mean direction and concentration parameters. An interesting property of this family is that vM distributions are fully defined by their first moment [6] \[\begin{align} m_1(\eta)&\triangleq \mathbb{E}(X)=\int_{x\in\mathcal{C}}{x \;\mathrm{vM}(x;\eta)} \mathrm{d}x \nonumber \\ & =A(\kappa) e^{j\mu}= A(|\eta|)\frac{\eta}{|\eta|}, \label{eq:vm95moment} \end{align}\tag{5}\] where \(A(\kappa)=\frac{I_1(\kappa)}{I_0(\kappa)}\), and \(I_i(.)\) denotes the \(i\)-th modified Bessel function of the first kind.

We now introduce a new iterative detection algorithm based on vM-BP [11] to perform message-passing-based MIMO detection on PSK-modulated symbols. The PSK modulation \(\Omega\) is a discrete subset of \(\mathcal{C}\); vM-BP consists of relaxing the domain of the optimization problem 2 from \(\Omega^K\) to \(\mathcal{C}^K\), to obtain simple parametric messages for the BP implementation.

3.2 Factor Graph Representation↩︎

Let us first consider the coherent case where the channel is perfectly known. The system model described in Sec. 2 with the factorization of 4 can be represented by the factor graph given in Fig. 1.

The graph contains two types of nodes: a) variable nodes (VNs), denoted as circular nodes, comprising information symbol nodes \(x^{(k)}\) and observation symbol nodes \(y_n\); b) function nodes (FNs), denoted as rectangular boxes. There are two types of function nodes, namely prior functions (denoted by \(f_{X^{(k)}}\)) and local likelihood functions attached to observations (denoted by \(f_n\)).

3.2.0.1 Priors

Starting with the prior functions, we assume that the true prior of information symbols is uniform on \(\Omega\), given by \(f_{X^{(k)}}(x) = \sum_{\alpha\in \Omega}\frac{1}{M} \delta(x-\alpha)\), where \(\delta\) is the Dirac delta function. Since vM-BP relies on the relaxation of the problem to the unit circle, we use instead a uniform prior on \(\mathcal{C}\), which can be expressed as a von Mises distribution with concentration parameter equal to zero, i.e. \(f_{X^{(k)}}(x) = \mathrm{vM}\left(x; 0\right)\).

3.2.0.2 Observations

The other type of function nodes consists of the local likelihood functions of the observations \({y_n}\) given \(\mathbf{x}\) and \(\mathbf{h}_n\), i.e. \[f_n\biggl(\mathbf{x}; y_n ,\mathbf{h}_n\biggr) = \frac{1}{\pi \sigma^2} \exp{ \biggl\{ \frac{-1}{\sigma^2} \left | y_n- \sum_{k} x^{(k)} h^{(k)}_n \right| ^2 \biggr\} }, \label{eq:likelihood95MIMO}\tag{6}\] where \(n \in \{1,\dots,N\}\), \(\mathbf{x} = [x^{(1)},\dots,x^{(K)}]\) and \(\mathbf{h}_n = [h^{(1)}_n,\dots,h^{(K)}_n]\). Given this representation, there are \(K\) symbols to be estimated based on \(N\) observations. Each information symbol node \(x^{(k)}\) is connected to \(N\) function nodes, while each function node \(f_n\) is connected to \(K\) information symbols. Thus, the total number of edges in the central part of the factor graph is equal to \(KN\).

3.3 vM-BP Messages and Beliefs↩︎

In this subsection we define an approximate belief propagation algorithm on the factor graph represented on Fig. 1, whereby all messages to and from an information symbol node are constrained to belong to the parametric family of von Mises distributions. This allows to replace a discrete probability distribution over \(\Omega\) by a continuous distribution over \(\mathcal{C}\) parameterized by a single complex parameter, as detailed in Section 3.1. We show in the sequel that this yields parametric expressions of the messages which also belong to the family of von Mises distributions, and furthermore admit simple parameter update expressions. Let us now consider the two types of messages that occur in the graph at iteration \(t\), and the resulting belief expressions.

3.3.1 Message from VN \(x^{(k)}\) to FN \(f_n\)↩︎

It is given by \[\lambda^{[t]}_{{x^{(k)} \rightarrow f_{n}} }(x^{(k)})= \prod_{{n'\neq n}} \lambda_{f_{{n'}} \rightarrow x^{(k)}}^{[t-1]}(x^{(k)}), \label{message95from95variable95node95to95function95node}\tag{7}\] where \(\lambda_{f_{{n'}} \rightarrow x^{(k)}}^{[t-1]}(x) \propto \textrm{vM}(\, x \,;\eta^{[t-1]}_{f_{n'} \rightarrow x^{(k)} })\). Therefore, the message \(\lambda^{[t]}_{x^{(k)} \rightarrow f_{n}}\) is a product of exponential distributions which can be identified as a von Mises distribution \(\lambda^{[t]}_{x^{(k)} \rightarrow f_{n}}(x) \propto \textrm{vM}(\,x\,;\eta^{[t]}_{x^{(k)} \rightarrow f_{n}})\), with natural parameter given as the sum of the natural parameters of the incoming messages at the previous iteration: \[\eta^{[t]}_{x^{(k)} \rightarrow f_{n}} = \sum_{{n'\neq n}} \eta_{f_{{n'}} \rightarrow x^{(k)}}^{[t-1]}. \label{eq:natural95parameter95vn95to95fn}\tag{8}\]

3.3.2 Message from FN \(f_n\) to VN \(x^{(k)}\)↩︎

This is given by \[\begin{align} \lambda_{f_{n} \rightarrow x^{(k)}}^{[t]}(x^{(k)}) = \int_{\mathbf{x}_{\sim x^{(k)}} \in \mathcal{C}^{K-1}} & f_n\biggl( \mathbf{x} ; y_n,\mathbf{h}_n\biggr) \\ & \prod_{k'\neq k} \lambda_{x^{({k'})} \rightarrow f_{n}}^{[t-1]}(x^{(k')}) d\mathbf{x}_{\sim x^{(k)}}, \\ \end{align} \label{eq:message95from95function95node95to95variable95node}\tag{9}\] where the integration is over the \(K-1\) information symbols from all users but the \(k^{\text{th}}\), denoted by the subscript \(\sim x^{(k)}\). Considering the definitions of the function node in 6 and of the messages \(\lambda^{[t]}_{x^{(k)} \rightarrow f_{n}}\), the integrand in 9 is a product of Gaussian and von Mises distributions, making the integral non-tractable. To overcome this issue, we make use of mean field approximation (see [9] and [10]) and replace 9 by \[\lambda_{f_{n} \rightarrow x^{(k)}}^{[t]} (x^{(k)}) = \exp{\left(\mathbb{E}_{\sim x^{(k)}}[\log(f_n(\mathbf{x}; y_n ,\mathbf{h}_n))]\right)}.\] Taking the logarithm of the message yields \[\begin{align} & \log \lambda_{f_{n} \rightarrow x^{(k)}}^{[t]}(x^{(k)}) = \\ & \frac{-1}{\sigma^2} \left| y_n- h_n^{(k)}x^{(k)} -\sum_{{k'}\neq k}m_1(\eta_{ x^{({k'})}\rightarrow f_n}^{[t-1]} )h^{({k'})}_n \right| ^2,\\ \end{align}\] where \(m_1(\eta_{ x^{({k'})}\rightarrow f_n}^{[t-1]} )\) is the first moment of the message from interfering variable node \(k'\) to factor node \(f_n\), as per 5 . After expanding and removing the constant terms with respect to \(x^{(k)}\) by using the fact \(\rvert x^{(k)}\rvert^2=1\), the message to \(x^{(k)}\) from function node \(f_n\) is given by \[\begin{align} & \lambda_{f_{n} \rightarrow x^{(k)}}^{[t]} (x^{(k)})\propto \\ & \exp{\biggl\{ \mathrm{Re} \left( \frac{2h_n^{(k)}}{\sigma^2}\left(y_n^* - \sum_{{k'}\neq k}m_1^*(\eta_{ x^{({k'})}\rightarrow f_n}^{[t-1]}) h^{*({k'})}_n \right) x^{(k)} \right) \biggr\}}.\\ \end{align} \label{eq:belief95fn95to95vn}\tag{10}\] This expression can be identified as a von Mises distribution with natural parameter \[\eta_{f_{n} \rightarrow x^{(k)}}^{[t]} = \frac{2(h_n^{(k)})^*}{\sigma^2}\left(y_n - \sum_{{k'}\neq k}m_1(\eta_{ x^{({k'})}\rightarrow f_n}^{[t-1]}) h^{({k'})}_n \right) \label{eq:natural95parameter95function95variable95final}.\tag{11}\]

3.3.3 Beliefs↩︎

The beliefs of information symbol \(x^{(k)}\) at iteration \(t\) can be computed according to \[b^{[t]}(x^{(k)}) \propto f_{X^{(k)}}(x^{(k)}) \prod_{n\in \{1,\dots,N\}} \lambda_{f_{{n}} \rightarrow x^{(k)}}^{[t]}(x^{(k)}). \label{belief95MIMO}\tag{12}\] Since every term in the product is a von Mises distribution, \(b^{[t]}(x^{(k)})\) is also von Mises, with natural parameter \[\hat{\eta}_{x^{(k)}}^{[t]} = \eta^{[0]}_{x^{(k)}} + \sum_{{n}} \eta_{f_{{n}} \rightarrow x^{(k)}}^{[t]}, \label{eta95MIMO}\tag{13}\] where \(\eta^{[0]}_{x^{(k)}}=0\) for the uniform prior.

3.4 vM-BP Algorithm Outline↩︎

At iteration 0 (initialization), there is no information but the prior, hence, \(\lambda^{[0]}_{x^{(k)} \rightarrow f_{n}}(x) = \text{vM}\left(\,x\,;0\right)\). Given this initial message, the first message from the function node (\(f_n\)) to the variable node (\(x^{(k)}\)) is \(\lambda_{f_{{n}} \rightarrow x^{(k)}}^{[0]}(x^{(k)}) = \text{vM}\left(\,x^{(k)}\,;2(h_n^{(k)})^*y_n/\sigma^2\right)\). The receiver then iteratively updates the messages in both directions using 8 and 11 for a fixed number of iterations \(N_{\mathrm{itr}}\), or until convergence. At the end of the algorithm, the information symbols are estimated from the beliefs 13 . The MAP hard decision rule for the estimated PSK symbol is given by \[\hat{x}^{(k)} = \arg\max_{x\in \Omega} \mathrm{vM}(x;\hat{\eta}_{x^{(k)}}^{[N_{\mathrm{itr}}]}), \;\forall k=1,\cdots,K. \label{estimation95x}\tag{14}\]

4 vM-BP detection for MIMO with imperfect channel knowledge↩︎

The discussion so far focused on estimation of information symbols \(x^{(k)}\) with perfectly known channel. Let us now consider a case where the receiver has an imperfect knowledge of the channel and tries to infer both the information symbols and the true channel state. Let \(\tilde{h}_n^{(k)}\) denote the imperfect channel knowledge available to the receiver. The prior of the channel variables captures the imperfect CSI through a complex Gaussian distribution given by \[\begin{align} f_{H_{n}^{(k)}}(h) & = \frac{1}{\pi \sigma_e^2} \exp{\biggl\{\frac{ -\lvert h -\mu_{h_n^{(k)}}^{[0]} \rvert ^2}{\sigma_e^2}\biggr\}}\\ & = \text{CN}\left(h;\mu_{h_n^{(k)}}^{[0]},\sigma_e^2\right), \end{align} \label{prior95pdf95h}\tag{15}\] with mean \(\mu_{h_n^{(k)}}^{[0]} = \tilde{h}_n^{(k)}\) and variance \(\sigma_e^2\).

Fig. 2 depicts the updated factor graph with an additional VN type for the channel state variables, denoted by \(h_n^{(k)}\), and an additional FN type, denoted by \(f_{H_n^{(k)}}\). In this representation, each channel state node is connected to a single FN, and observation FNs \(f_n\) are connected to \(2K\) VNs. In this factor graph, in Fig. 2, there are four types of messages that can occur: messages to and from \(x^{(k)}\) nodes; and messages to and from \(h_n^{(k)}\) nodes. The messages to and from a channel node are constrained to be complex normal (CN) distributions; and the messages to and from an information symbol node are constrained to be vM distributions. The message from VN \(x^{(k)}\) to FN \(f_n\) remains the same as in Sec. 3.3.1. The rest of the messages are discussed below.

4.1 Message from VN \(h_n^{(k)}\) to FN \(f_n\)↩︎

This is given by \[\lambda^{[t]}_{h_n^{(k)} \rightarrow f_{n}}(h_n^{(k)}) = \prod_{n'\neq n} \lambda_{f_{{n'}} \rightarrow h_n^{(k)}}^{[t-1]}(h_n^{(k)}), \label{message95h95to95fn95imperfec95def}\tag{16}\] where \(\lambda_{f_{{n'}} \rightarrow h_n^{(k)}}^{[t-1]} (h_n^{(k)}) \propto \text{CN}\left(h_n^{(k)}; \mu^{[t-1]}_{f_{n'} \rightarrow h_n^{(k)} },\sigma^{2[t-1]}_{f_{n'} \rightarrow h_n^{(k)} } \right)\) and \(n'\) takes all values in \(\{1,\cdots,N\}\) that are connected to the \(h_n^{(k)}\) node except \(n\). Since every term in 16 is a CN distribution, the resultant is also a CN distribution given as \(\lambda^{[t]}_{h_n^{(k)} \rightarrow f_{n}}(h_n^{(k)}) = \text{CN}\left(h_n^{(k)}; \mu^{[t]}_{h_n^{(k)} \rightarrow f_{n}},\sigma^{2[t]}_{h_n^{(k)} \rightarrow f_{n}} \right)\), with parameters given as \[\begin{align} & \sigma^{2[t]}_{h_n^{(k)} \rightarrow f_{n}} = \left(\sum_{n'}\frac{1}{\sigma^{2[t-1]}_{f_{n'} \rightarrow h_n^{(k)}}}\right)^{-1},\\ & \mu^{[t]}_{h_n^{(k)} \rightarrow f_{n}} = \sigma^{2[t]}_{h_n^{(k)} \rightarrow f_{n}} \sum_{n'}\left( \frac{\mu^{[t-1]}_{f_{n'} \rightarrow h_n^{(k)}}}{\sigma^{2[t-1]}_{f_{n'} \rightarrow h_n^{(k)}}} \right).\\ \end{align} \label{parameters95h95to95f95imperfec95chan95final}\tag{17}\]

4.2 Message from FN \(f_n\) to VN \(x^{(k)}\)↩︎

This is given by \[\begin{align} & \lambda_{f_{n} \rightarrow x^{(k)}}^{[t]}(x^{(k)}) = \int_{\mathbf{x}_{\sim x^{(k)}} \in \mathcal{C}^{K-1}, \mathbf{h}_n \in \mathbb{C}^K} f_n\biggl( \mathbf{x},\mathbf{h}_n ; y_n\biggr)\\ & \prod_{k'\neq k} \lambda_{x^{({k'})} \rightarrow f_{n}}^{[t-1]}(x^{(k')}) \prod_{k} \lambda_{h_n^{({k'})} \rightarrow f_{n}}^{[t-1]}(h_n^{(k')}) d\mathbf{x}_{\sim x^{(k)}} d\mathbf{h}_n, \\ \end{align} \label{eq:message95from95function95node95x95to95variable95node95imperfect}\tag{18}\] where the integration is over \(K-1\) information symbols and \(K\) channel nodes. Using the same technique as in Sec. 3.3.2 to solve this integral, the message becomes \(\lambda_{f_{n} \rightarrow x^{(k)}}^{[t]}(x^{(k)}) = \text{vM}\left(x^{(k)};\eta_{f_{n} \rightarrow x^{(k)}}^{[t]}\right)\), where \[\begin{align} & \eta_{f_{n} \rightarrow x^{(k)}}^{[t]} = \\ & \frac{2(\mu_{h_n^{({k})} \rightarrow f_{n}}^{[t-1]})^*}{\sigma^2}\left(y_n - \sum_{{k'}\neq k}m_1(\eta_{ x^{({k'})}\rightarrow f_n}^{[t-1]}) \mu_{h_n^{({k'})} \rightarrow f_{n}}^{[t-1]} \right).\\ \end{align}\]

4.3 Message from FN \(f_n\) to VN \(h_n^{(k)}\)↩︎

This is given by \[\begin{align} & \lambda_{f_{n} \rightarrow h_n^{(k)}}^{[t]}(h_n^{(k)}) = \int_{{\mathbf{h}_n}_{\sim h_n^{(k)}} \in \mathbb{C}^{K-1}, \mathbf{x} \in \mathcal{C}^K} f_n\biggl( \mathbf{h}_n,\mathbf{x} ; y_n\biggr)\\ & \prod_{k'\neq k} \lambda_{h_n^{({k'})} \rightarrow f_{n}}^{[t-1]}(h_n^{(k')}) \prod_{k} \lambda_{x^{({k'})} \rightarrow f_{n}}^{[t-1]}(x^{(k')}) d\mathbf{x} d{\mathbf{h}_n}_{\sim h_n^{(k)}}, \\ \end{align} \label{eq:message95from95function95node95x95to95variable95node}\tag{19}\] where the integration is over \(K-1\) channel nodes and \(K\) information symbol nodes. Using the same technique as in Sec. 3.3.1 to solve this integral, the message becomes \[\lambda_{f_{n} \rightarrow h_n^{(k)}}^{[t]}(h_n^{(k)}) = \text{CN} \left(h_n^{(k)};\mu_{f_{n} \rightarrow h_n^{(k)}}^{[t]}, \sigma_{f_{n} \rightarrow h_n^{(k)}}^{2[t]}\right),\] where the parameters are written as \[\begin{align} &\mu_{f_{n} \rightarrow h_n^{(k)}}^{[t]} \\ & = m_1(\eta_{ x^{({k})}\rightarrow f_n}^{[t-1]}) \left( y_n^* - \sum_{k'\neq k}m_1^*(\eta_{ x^{({k'})}\rightarrow f_n}^{[t-1]}) \mu_{ h_n^{({k'})}\rightarrow f_n}^{*[t-1]}\right),\\ & \sigma_{f_{n} \rightarrow h_n^{(k)}}^{2[t]} = \sigma^2. \end{align}\]

4.4 Beliefs and Estimation↩︎

The beliefs for VN \(x^{(k)}\) are calculated using 13 and the corresponding estimation of symbols from 14 . Similarly, the belief of VN \(h_n^{(k)}\) can be computed using \[b^{[t]}(h_n^{(k)}) \propto f_{H_n^{(k)}}(h_n^{(k)}) \prod_{n'} \lambda_{f_{{n}} \rightarrow h_n^{(k)}}^{[t]}(h_n^{(k)}), \label{belief95chan}\tag{20}\] where \(n'\) takes all values in \(\{1,\cdots,N\}\) that are connected to the \(h_n^{(k)}\) node. Since all terms are CN distribution, the belief of VN \(h_n^{(k)}\) is also a CN distribution, written as \(b^{[t]}(h_n^{(k)}) = \text{CN} \left( h_n^{(k)}; \mu_{h_n^{(k)}}^{[t]},\sigma_{h_n^{(k)}}^{2[t]}\right)\), where the parameters are given by \[\begin{align} & \sigma^{2[t]}_{h_n^{(k)}} = \left(\frac{1}{\sigma_h^{2[0]}} + \sum_{n'} \frac{1}{\sigma^{2[t]}_{f_{n'} \rightarrow h_n^{(k)}}}\right)^{-1},\\ & \mu^{[t]}_{h_n^{(k)}} = \sigma^{2[t]}_{h_n^{(k)}} \left(\frac{\mu^{[0]}_{h_n^{(k)}}}{\sigma^{2[0]}_{h}} + \sum_{n'} \frac{\mu^{[t]}_{f_{n'} \rightarrow h_n^{(k)}}}{\sigma^{2[t]}_{f_{n'} \rightarrow h_n^{(k)}}} \right).\\ \end{align}\] The channel state is estimated from the belief and the estimated symbol is given by \[\hat{h}_n^{(k)} = \mu^{[N_{itr}]}_{h_n^{(k)}}.\]

5 Performance Analysis↩︎

5.1 Complexity for Coherent Case↩︎

Let us first evaluate the complexity of the proposed algorithm in terms of floating point operations (FLOPs) per iteration, following the per-operation complexity from [14]. In vM-BP, instead of tracking a mean and a variance as usual in approximate inference algorithms, message updates can be derived using only the complex natural parameter \(\eta\).

The naive complexity number obtained by multiplying the respective complexities of computing 8 and 11 by the number of their occurrences at each iteration of the algorithm (\(KN\) times, one message per edge) can be reduced by noticing that the sums of all-but-one terms in 8 and 11 which must be computed for all \(k\) and \(n\) have many terms in common. Computing the total sum once per iteration, either per observation or per user and removing the contribution from the variable of interest to compute each message leads to a complexity of \(33KN-10N-2K+KNC_A\) FLOPS, where \(C_A\) is the cost of computing \(A(\cdot)\). The efficient evaluation of \(A(\cdot)\) is a rich research topic [15] which we do not address in this article; note however that it is a well-behaved scalar function of a scalar parameter, which can be approximated efficiently. Hence we keep \(C_{A}\) undefined for the purpose of complexity evaluation. Note that the complexity of vM-BP is independent of the modulation order \(M\) since the relaxation to the unit circle obliterates the need to consider the discrete modulation \(\Omega\) during the iterations.

In Table 1, the complexity of the proposed vM-BP is compared to three other iterative algorithms: EP (for which we also assume the use of the aforementioned complexity reduction approach) and generalized-AMP (G-AMP) [5], which are both based on the same factor graph representation as vM-BP, as well as VEP [13], in which all receive antennas are handled by a single factor node. VEP generally performs better at the cost of increased computational cost because each iteration involves a matrix inversion.

In Fig. 3, the number of FLOPs is plotted against the number of users \(K\) with the number of antennas set as \(N = 4K\) (lightly loaded system), for different values of the PSK order \(M\). Different iterative algorithms like VEP, EP, GAMP and vM-BP are compared. It can be seen that VEP is more complex compared to the rest of the decoders.

None

5.2 Error rate analysis↩︎

In this section, we evaluate the performance of the proposed vM-BP for coherent MIMO detection in terms of symbol error rate (SER) for different system parameters. We normalize each channel vector in 1 such that \(||\mathbf{h}^{(k)}||=1\) for all \(k\), and first assume exact knowledge of the channel state. The signal to noise ratio (SNR) is defined as \(1/\sigma^2\). For simulations, a system with \(K=16\) users and \(N=64\) receive antennas is considered. The proposed vM-BP algorithm is compared with L-MMSE, EP, AMP and VEP. In this scenario with \(K\ll N\), we consistently observed that damping had almost no impact, hence it is not considered in this section.

Let us first focus on the convergence behavior. Fig. 4 depicts the SER versus number of iterations for SNR = 16 dB and \(M=8\) for vM-BP, G-AMP and EP, for Rayleigh fading. All algorithms require about \(5\) iterations to converge. EP has the lowest plateau with a minimal SER equal to \(7.10^{-4}\), G-AMP stops at \(9.10^{-4}\) and vM-BP at \(10^{-3}\). In the sequel, the number of iterations is set to \(10\) to ensure convergence.

None

Performance in terms of SER versus SNR is evaluated over two different channel models: a) Gaussian i.i.d. channel and b) Correlated channel modeled using angle of arrivals (AoA).

5.2.1 Gaussian i.i.d. Channels↩︎

Fig. 5 shows the SER against SNR for different \(M\) assuming that the receiver has perfect knowledge of the true Gaussian i.i.d. (Rayleigh fading) channel. The performance of VEP is not depicted because it performs identically to G-AMP in this scenario.

None

The three iterative algorithms outperform L-MMSE detection for all considered modulation orders (\(M=8\), 32, 64). While AMP outperforms vM-BP at low SNR, which is coherent with Fig. 4, the gap becomes smaller as the SNR increases. Scalar EP consistently outperforms vM-BP by \(0.5\) dB.

5.2.2 Correlated channels↩︎

In this section, the channels follow a geometric model assuming line-of-sight propagation and a uniform linear array with half-wavelength spacing at the receiver side, hence the channel of each user consists in the array steering vector which is fully determined by the angle of arrival (AoA). We consider three different scenarios, which differ in how the AoAs are chosen: (i) high correlation, with AoAs for user \(k\) given as \(k°\) from the broadside direction, where the channel matrix \(\mathbf{H}\) has condition number as \(4.6 \cdot 10^4\), (ii) medium correlation with AoAs linearly spaced between \([-10°, 10°]\), with channel condition number \(201\), and (iii) low correlation, with AoAs linearly spaced in \([-60°, 60°]\) for which the channel condition number is \(1.2\).

None

The performance curves depicted in Fig. 6 show that, for all three cases, vM-BP performs better than L-MMSE. As expected, increased correlation (and channel condition number) is detrimental to all detection algorithms. All detectors exhibit the same performance for the low correlation scenario. This is expected since in LoS conditions with hardly any correlation with the other users the detection is easy. The performance of G-AMP is not depicted, it is essentially identical to that of EP in this scenario. In the high correlation scenario, only VEP is able to detect somewhat reliably. The results for the medium correlation are different. The error floor for vM-BP, at \(8.10^{-3}\), is lower than the one for scalar EP, at \(3.10^{-2}\). Also, although VEP does not suffer from an error floor, it is outperformed by vM-BP from 13 dB to 23 dB. No intuitive explanation is available for the discrepancy between the performance rankings of the various algorithms over high and medium correlated channels.

5.2.3 Overloaded System with Gaussian i.i.d. Channel↩︎

We consider an overloaded system with \(K = 120\) users and \(N = 96\) receive antennas, for PSK of order \(M = 8\), with Gaussian i.i.d. channel. The SER versus number of iterations for SNR = 16 dB is plotted in Fig. 7 for different damping factors. This damping is introduced in the messages from FN to VN and an improvement in performance is noted. All algorithms seem to converge at around 25 iterations. Hence we fix the number of iterations to 25 for all algorithms in the next simulation, depicted in Fig. 8 where SER is plotted against SNR.

None

For Fig. 8, the performance is compared with L-MMSE, G-AMP, VEP and EP. As expected, the iterative algorithms significantly outperforms L-MMSE, which as a linear method, is not able to resolve the over-determined system of equations even at high SNR. From Fig. 8 it is evident that the iterative algorithms do not perform well with no damping and an improvement in performance is noticed when damping is introduced.

None

5.3 Imperfect Channel Knowledge↩︎

In this subsection, we analyze the performance of the algorithm from Sec. 4 when only imperfect CSI is available. The imperfect channel is modeled as \[\tilde{h}_n^{(k)} = h_n^{(k)} + e, \label{imperfect95chan95mod}\tag{21}\] where the true channel \(h_n^{(k)}\) is Gaussian i.i.d. (with no normalization) and \(e\) is a complex channel estimation error with zero mean and variance \(\sigma_e^2\). We consider three scenarios: (i) a mismatched detector, where the receiver assumes that the true CSI is equal to \(\tilde{h}_n^{(k)}\) and works on the factor graph from Fig. 1 to estimate only the information symbols; (ii) a CSI-imperfection-aware detector where the receiver knows that \(\tilde{h}_n^{(k)}\) is imperfect CSI and is aware of the estimation error variance \(\sigma_e^2\), and tries to obtain an improved estimate of CSI along with information symbols, using the factor graph from Fig. 2; (iii) a block-fading scenario with (known) blocklength \(L\), whereby each user transmits \(L\) symbols while the channel remains unchanged; the receiver tries to estimate both the symbols and channel. This is an extension of scenario (ii) where the graph is extended by adding information symbol nodes. SER performance curves are depicted in Fig. 9. For the mismatched decoder case we plot L-MMSE performance, while for the imperfection-aware decoder we plot vM-BP. For \(L = 1\), i.e., each user transmitting only one symbol, damping did not improve the performance. However, for \(L > 1\), damping showed significant improvement in performance, and we fix the damping coefficient to 0.6 for messages from FN to VN (both channel and information symbol nodes). The number of iterations for this scenario is set to 20.

None

It can be noted that, the imperfection-aware detector improves the performance. As the number of symbols transmitted from each user increase, SER decreases, as the channel estimates are improved since the degree of channel node increases. As \(L\) increases, the SER performance approaches the performance achievable under perfect CSI, even when provided with a highly imperfect prior (high \(\sigma_e^2\)).

6 Conclusion↩︎

This paper proposed a BP detector for mMIMO with PSK modulations based on von Mises distributions by leveraging directional statistics and relaxing the discrete constellation on the unit circle to a continuous vM distribution. This results in an algorithm that is fully continuous and agnostic to the modulation order. The numerical results and complexity of the algorithm are discussed in detail and compared to state of the art iterative detectors, such as EP, G-AMP. Extensions of the proposed algorithm to the presence of imperfect channel are also presented, making the proposed algorithm suitable for practical mMIMO systems.

References↩︎

---

1. This work has been partially supported by the ANR under the “France 2030” program under grant NF-PERSEUS (ANR-22-PEFT-0004), and under grant WARM-M2M (ANR-24-CE25-2514).↩︎