Expectation Propagation Detector for Extra-Large Scale Massive MIMO

TL;DR

This paper introduces a subarray-based expectation propagation detector for extra-large-scale massive MIMO systems, balancing performance and practicality by reducing complexity and enabling distributed processing.

Contribution

It extends the expectation propagation detector to subarray architectures, providing convergence analysis, complexity reduction strategies, and demonstrating improved performance over existing methods.

Findings

The proposed detector converges within few iterations.

It outperforms existing detectors in simulations.

The method reduces complexity and data exchange overhead.

Abstract

The order-of-magnitude increase in the dimension of antenna arrays, which forms extra-large-scale massive multiple-input-multiple-output (MIMO) systems, enables substantial improvement in spectral efficiency, energy efficiency, and spatial resolution. However, practical challenges, such as excessive computational complexity and excess of baseband data to be transferred and processed, prohibit the use of centralized processing. A promising solution is to distribute baseband data from disjoint subsets of antennas into parallel processing procedures coordinated by a central processing unit. This solution is called subarray-based architecture. In this work, we extend the application of expectation propagation (EP) principle, which effectively balances performance and practical feasibility in conventional centralized MIMO detector design, to fit the subarray-based architecture. Analytical results confirm the convergence of the proposed iterative procedure and that the proposed detector asymptotically approximates Bayesian optimal performance under certain conditions. The proposed subarray-based EP detector is reduced to centralized EP detector when only one subarray exists. In addition, we propose additional strategies for further reducing the complexity and overhead of the information exchange between parallel subarrays and the central processing unit to facilitate the practical implementation of the proposed detector. Simulation results demonstrate that the proposed detector achieves numerical stability within few iterations and outperforms its counterparts.

Tables

Table 1. TABLE I: Complexity Comparison

Algorithm		Each Parallel LPM		CPM
		Full ${𝐇}_c$	Trimmed ${\widetilde{𝐇}}_c$
Algorithm 1	Mult.	$8N_{c}TK\left(K+1\right)+$	$8N_{c}T{K}_c\left(K_c+1\right)+$	$CT\left(1+2K\right)+$
		$6TK\left(K+2\right)$	$6T{K}_c\left(K_c+2\right)$	$K\left(T-1\right)\left(7\left|𝒳\right|+2\right)$
	Exp.	0	0	$K\left|𝒳\right|\left(T-1\right)$
	Trans.	$T\left(2K+1\right)$	$T\left(2K_c+1\right)$	$\left(2K+1\right)\left(T-1\right)$
One feedforward	Mult.	$TK\left[8N_c\left(K+1\right)+2\left(4K+3\right)\right]$	$T{K}_c\left[8N_c\left(K_c+1\right)+2\left(4K_c+3\right)\right]$	$C\left(1+2K\right)$
		$+K\left(T-1\right)\left(7\left|𝒳\right|+6\right)$	$+K_c\left(T-1\right)\left(7\left|𝒳\right|+6\right)$
	Exp.	$KT\left|𝒳\right|$	$K_{c}T\left|𝒳\right|$	0
	Trans.	$2K+1$	$2K_c+1$	0
Centralized EP	Mult.	$TK\left[8N\left(K+1\right)+2\left(4K+3\right)\right]+K\left(T-1\right)\left(7\left|𝒳\right|+6\right)$
	Exp.	$K\left|𝒳\right|\left(T-1\right)$

Equations

$$\mathbf{y}=\mathbf{Hx}+\mathbf{n},$$

$${\mathbf{y}}_c={\mathbf{H}}_c\mathbf{x}+{\mathbf{n}}_c,$$

$${\mathbf{y}}_c={\tilde{\mathbf{H}}}_c{\mathbf{x}}_c+{\mathbf{n}}_c,$$

$$p(\mathbf{x}|\mathbf{y},\mathbf{H})\propto \underset{\propto p(\mathbf{y}|\mathbf{x},\mathbf{H})}{\underbrace{\exp \left(-\frac{\Vert \mathbf{y}-\mathbf{Hx}{\Vert }^2}{{\sigma }^2}\right)}}\underset{=p(\mathbf{x})}{\underbrace{\prod _{k=1}^{K}p(x_k)}},$$

$$p(\mathbf{x}|\mathbf{y},\mathbf{H})\propto f_0(\mathbf{x})\prod _{c\in \mathcal{C}}{f}_c(\mathbf{x}),$$

$${\tau }_c={\omega }_0-{\eta }_c$$

$${\mathbf{\gamma }}_c={\tau }_c^{-1}\left(\frac{{\omega }_0{\hat{\mathbf{x}}}_0-{\eta }_c\text{buildrel}\text{lower}3.0pt\to }{{\mathbf{x}}_c}\right)$$

$${\mathbf{\Sigma }}_c={\left({\sigma }^{-2}{\mathbf{H}}_c^H{\mathbf{H}}_c+{\tau }_c{\mathbf{I}}_K\right)}^{-1}$$

$${\hat{\mathbf{x}}}_c={\mathbf{\Sigma }}_c\left({\sigma }^{-2}{\mathbf{H}}_c^H{\mathbf{y}}_c+{\tau }_c{\mathbf{\gamma }}_c\right)$$

$${\omega }_c={\left(K^{-1}\mathsf{tr}\left({\mathbf{\Sigma }}_c\right)\right)}^{-1}$$

$${\eta }_c={\omega }_c-{\tau }_c$$

$$\frac{\text{buildrel}\text{lower}3.0pt\to }{{\mathbf{x}}_c={\eta }_c^{-1}\left({\omega }_c{\hat{\mathbf{x}}}_c-{\tau }_c{\mathbf{\gamma }}_c\right)}$$

$${\tau }_0=\sum _{c=1}^C{\eta }_c$$

$$\frac{{\mathbf{\gamma }}_0={\tau }_0^{-1}\sum _{c=1}^C{\eta }_c\text{buildrel}\text{lower}3.0pt\to }{{\mathbf{x}}_c}$$

$${\hat{x}}_{0,k}=\mathsf{E}\left\{x_k|{\gamma }_{0,k},{\tau }_0\right\},\text{for}k\in \mathcal{K}$$

$$v_{0,k}=\mathsf{Var}\left\{x_k|{\gamma }_{0,k},{\tau }_0\right\},\text{for}k\in \mathcal{K}$$

$${\omega }_0={\left(K^{-1}\sum _{k=1}^K{v}_{0,k}\right)}^{-1}$$

$${\mathbf{\gamma }}_0=\mathbf{x}+{\mathbf{n}}_0,$$

$$\mathrm{P}(x_k|{\gamma }_{0,k},{\tau }_0)=\frac{\mathcal{C}\mathcal{N}({\gamma }_{0,k};x_k,{\tau }_0^{-1})\mathrm{P}(x_k)}{\sum _{x\in \mathcal{X}}\mathcal{C}\mathcal{N}({\gamma }_{0,k};x_k,{\tau }_0^{-1})\mathrm{P}(x)},$$

$${\hat{x}}_{0,k}$$

$$v_{0,k}$$

$$\frac{{\mathbf{\gamma }}_0=\frac{1}{{\tau }_0}{\sum }_{c=1}^C\text{buildrel}\text{lower}3.0pt\to }{{\mathbf{x}}_c,{\mathbf{\gamma }}_c=\frac{1}{{\tau }_c}\left(\frac{{\omega }_0{\hat{\mathbf{x}}}_0-\text{buildrel}\text{lower}3.0pt\to }{{\mathbf{x}}_c}\right).}$$

$${\rho }^t={\left(\sum _{c=1}^C\left(\frac{1}{{\mathsf{mse}}_c\left({\nu }_c^t\right)}-\frac{1}{{\nu }_c^t}\right)\right)}^{-1},$$

$${\nu }_c^{t+1}={\left(\frac{1}{{\mathsf{mse}}_0\left({\rho }^t\right)}-\left(\frac{1}{{\mathsf{mse}}_c\left({\nu }_c^t\right)}-\frac{1}{{\nu }_c^t}\right)\right)}^{-1}\mathrm{for}c\in \mathcal{C},$$

$${\mathsf{mse}}_c\left({\nu }_c^t\right)=\frac{1}{K}\sum _{i=1}^K\frac{{\sigma }^2{\nu }_c^t}{{\lambda }_{i,c}{\nu }_c^t+{\sigma }^2},$$

$${\mathsf{mse}}_0\left({\rho }^t\right)={\mathsf{E}}_r\left\{{\mathsf{E}}_{x|r}\left[{\left|x-{\mathsf{E}}_{x|r}\left[x|r\right]\right|}^2\right]\right\}.$$

$${\omega }_0={\omega }_1=\cdots ={\omega }_C=\omega :=\frac{1}{C}\left({\tau }_0+\sum _{c=1}^C{\tau }_c\right),$$

$${\hat{\mathbf{x}}}_0={\hat{\mathbf{x}}}_1=\cdots ={\hat{\mathbf{x}}}_C=\hat{\mathbf{x}}:=\frac{1}{C}\cdot \frac{{\tau }_0{\mathbf{\gamma }}_0+{\sum }_{c=1}^C{\tau }_c{\mathbf{\gamma }}_c}{{\tau }_0+{\sum }_{c=1}^C{\tau }_c}.$$

$$b_0(\mathbf{x}):=\frac{1}{Z_0\left({\mathbf{\gamma }}_0\right)}p(\mathbf{x})\exp \left[-{\tau }_0\Vert \mathbf{x}-{\mathbf{\gamma }}_0{\Vert }^2\right],$$

$$b_c(\mathbf{x}):=\frac{1}{Z_c\left({\mathbf{\gamma }}_c\right)}\exp \left[-\frac{1}{{\sigma }^2}\Vert {\mathbf{y}}_c-{\mathbf{H}}_c\mathbf{x}{\Vert }^2-{\tau }_c\Vert \mathbf{x}-{\mathbf{\gamma }}_c{\Vert }^2\right],$$

$$q(\mathbf{x}):=\frac{1}{Z_q\left(\hat{\mathbf{x}}\right)}\exp \left[-\omega \Vert \mathbf{x}-\hat{\mathbf{x}}{\Vert }^2\right],$$

$$\underset{b_0,b_1,\dots ,b_C}{\min }\underset{q}{\max }J\left(b_0,b_1,\dots ,b_C,q\right)$$

$$:=\sum _{c=0}^C\mathsf{KL}\left(b_c\Vert e^{\log f_c(\mathbf{x})}\right)+CH(q)$$

$$\text{s.t.}{\mathsf{E}}_{b_c}(\mathbf{x})={\mathsf{E}}_q(\mathbf{x}),$$

$$\frac{1}{K}\sum _{k\in \mathcal{K}}{\mathsf{E}}_{b_c}(|x_k{|}^2)=\frac{1}{K}\sum _{k\in \mathcal{K}}{\mathsf{E}}_q(|x_k{|}^2),$$

$$\hat{\mathbf{x}}={\mathsf{E}}_{b_0}(\mathbf{x})={\mathsf{E}}_{b_1}(\mathbf{x})=\cdots ={\mathsf{E}}_{b_C}(\mathbf{x})={\mathsf{E}}_q(\mathbf{x}),$$

$$\begin{array}{rl}{\omega }^{-1}& =\frac{1}{K}\sum _{k\in \mathcal{K}}{\mathsf{E}}_{b_0}(|x_k{|}^2)=\frac{1}{K}\sum _{k\in \mathcal{K}}{\mathsf{E}}_{b_1}(|x_k{|}^2)=\cdots \\ & =\frac{1}{K}\sum _{k\in \mathcal{K}}{\mathsf{E}}_{b_C}(|x_k{|}^2)=\frac{1}{K}\sum _{k\in \mathcal{K}}{\mathsf{E}}_q(|x_k{|}^2).\end{array}$$

$$p\left({\mathbf{\gamma }}_0|\mathbf{x}\right)=\frac{p\left({\mathbf{\gamma }}_0\right)}{Z_0\left({\mathbf{\gamma }}_0\right)}\exp \left[{\tau }_0\Vert \mathbf{x}-{\mathbf{\gamma }}_0{\Vert }^2\right],$$

Full text

I Introduction

Massive multiple-input multiple-output (MIMO) is widely believed to be one of the key techniques for fifth-generation wireless systems [1]. This technique was proposed to equip a base station (BS) with a large-scale antenna array (in the order of hundreds or thousands) and simultaneously serve a relatively small number of mobile users in the same time-frequency resource [2]. Many advantages, including substantial improvement in spectral efficiency, energy efficiency, and spatial resolution, can be achieved by increasing the antenna array dimension [3]. This approach is expected to facilitate ultra high data rate and system throughput. These advantages can be obtained by using a simple linear transceiver [4]. Massive MIMO communication systems have recently been implemented in the academia and various industries, and many prototyping systems have been constructed [5, 6], which implies the potential commercial success of these communication systems in the future.

Current cellular networks commonly deploy compact and co-located antenna array with small antenna separation, such as in the order of the wavelength, on the top of a tower or roof. However, the conventional compact antenna deployment encounters numerous practical challenges when the array dimension is increased to the order of thousands and more, such as the size and weight of the array and wind load. Moreover, to achieve great spatial resolution, it is desirable to distribute antennas over a substantially large area [7]. These factors advocate embracing new antenna deployment strategies for arrays with extremely large dimensions. One potential approach is to integrate the antenna array into large structures, such as along the walls of tall buildings, airports, or large shopping malls, or along the structure of a stadium [8]. Another approach is to distribute disjoint subsets of antennas of the entire array over a large geographical area coordinated by a central processing unit, thus yielding distributed antenna systems [9]. We refer to communication paradigms that utilize antenna arrays with the order-of-magnitude increase in their dimensions as extra-large-scale massive MIMO systems.

In extra-large-scale massive MIMO systems, the conventional processing architecture, in which signals received by all antennas and the full channel matrix are involved in the centralized processing, faces numerous challenges, especially in crowded scenarios. An evident challenge is the excessive computation complexity due to the large number of antennas and users. Even simple linear transceivers involve complex matrix operations to high-dimensional matrices. Another prominent challenge is the need to transfer excessively large amounts of baseband data received by extra-large-scale antenna array to the baseband processing unit [10]. A promising solution is to divide the entire antenna array into a few disjoint units, referred to as subarrays. Each subarray is associated with an individual processing unit that accesses only its local signals [11] and performs parallel processing with reduced complexity to produce coarse estimate of the transmitted signal. A central processing unit is responsible for producing refined estimate by performing a certain combination to the outputs of parallel subarrays. This subarray-based architecture enables the use of parallel computation supporting components (e.g., graphics processing unit) and relaxes the bandwidth requirement of the interface circuits by deploying parallel interconnections between subarrays and their local processing units. In practical systems, subarrays may correspond to separate hardware entities or software-defined logical interconnections between different portions of a signal array and their flexibly assigned processing resource [12]. The number and size of subarrays are fixed in the former, whereas they are adjustable in the latter.

Spatial non-stationary channel property generally occurs on large-scale antenna arrays [11]. The term “non-stationary” means that the different parts of the array may observe the same channel paths with varying power or distinct channel paths [13, 14]. Particularly, energy received from each user by different portions of the array varies. In this case, majority of energy received a specific user concentrates on small portions of the entire array, thus, inherent sparsity is introduced in the spatial dimension of the channel matrix [12]. Meanwhile,“spatial stationary” means that the entire array receives approximately the same amount of energy from each user, and it usually features MIMO channels where a moderate number of (several tens of) antennas are compactly deployed. When spatial non-stationarity occurs, this subarray-based processing architecture becomes favorable. First, each subarray only needs to detect signals transmitted by users that have sufficiently large received energy on that subarray, thereby enabling complexity reduction. Moreover, spatial non-stationary channel properties complicates the modeling of the channel on the entire array, which harms accurate channel estimation. Nevertheless, simplification of channel acquisition can be realized by decomposing the entire array into subarrays whose channels are approximated as stationary [15], which inspires detector development based on subarrays.

This study focuses on the development of a subarray-based multiuser detector for extra-large-scale massive MIMO systems. A few preliminary studies have been conducted on this topic. The performance of the linear receiver for such a scenario was evaluated in [16]. The subsequent work [12] adopted the successive interference cancellation (SIC) principle, which detects a signal from a specific user from a subarray with favorable interference conditions and then removes its contributions from the other subarrays to achieve performance improvement. The main drawback of this work is that detection of a certain user is performed from one subarray only, therefore, the contributions from other subarrays to the signal detection of this user cannot be combined. Moreover, detection of different users is performed in a serial manner, which does not fully exploit the parallel computation potential of subarray-based architecture. These drawbacks motivate the use of advanced techniques to pursue better design. The researchers in [17] proposed a MIMO detector based on expectation propagation (EP) and conventional centralized processing architecture and reported that the proposed detector outperforms many classic and state-of-the-art solutions with moderate computational complexity. An integrated circuit design for the detector proposed in [17] was presented in [18], which showed that EP-based design effectively balances performance and practical feasibility. These facts inspire us to extend the application of EP to subarray-based processing architecture. The proposed detector is constructed based on the general form of the EP principle [19], and several analytical results are derived. Moreover, additional strategies are proposed to reduce complexity and information exchange overhead, and improve parallelism, to facilitate the practical implementation of the proposed design. We derive the following meaningful findings and insights, which could provide guidance to further studies.

• •

Algorithm derivation: The proposed detector is derived via message updating and passing based on EP principle [20, 19, 21] in a factor graph that describes the subarray-based architecture. This derivation procedure reveals the underlying logic and rationale of the proposed algorithm and makes it more than a purely heuristic one. Many extant methods, such as [17, 22, 23, 24], can be derived using the same derivation framework, which implies that the EP principle will continue to help in developing effective algorithms by describing a particular application scenario with a factor graph on which feasible computation of messages is allowed.

• •

Theoretical analyses: Convergence of the proposed iterative procedure is intuitively verified through evolution analysis and is justified by simulation results. We prove that the asymptotic behavior of the proposed subarray-based EP detector matches the replica prediction of Bayesian optimal performance when the channel matrix has independent and identically distributed (i.i.d.) entries. The simulation results reveal that even with a non-i.i.d. channel matrix, the performance under different sizes and numbers of subarrays is similar to that under centralized EP detector [17]. This finding enables complexity minimization by including a small number of antennas in each subarray without substantially performance degradation.

• •

Engineering efforts: We also propose several strategies to facilitate efficient implementation of the proposed detector. Specifically, we demonstrate that the sparsity in the local channel matrix of each subarray derived from the non-stationarity allows for additional complexity reduction. Furthermore, a hierarchical implementation architecture is proposed to decompose the operations in each subarray into several parallel computations, thus improving parallelism. A modification of the proposed detector is presented. In such modification, only one feedforward from each subarray to the central processing unit is required in the entire procedure, which remarkably decreases the information exchange overhead and latency.

Notations: In this paper, we use lowercase and uppercase boldface letters to represent vectors and matrices, respectively. For vector $\mathbf{a}$, $a_{j}$ denotes the $j$-th entry of $\mathbf{a}$; $\Re{\mathbf{a}}$ and $\Im{\mathbf{a}}$ denote its real and imaginary parts respectively; and $[\mathbf{a}]_{{\mathcal{A}}}$ denotes the vector comprised by the elements of $\mathbf{a}$ indexed by ${\mathcal{A}}$. For matrix $\mathbf{A}$, we indicate its conjugate transpose, trace and inverse by $\mathbf{A}^{H}$, ${\sf tr}(\mathbf{A})$ and $\mathbf{A}^{-1}$ respectively. The operator $\odot$ denotes the element-wise (Hadamard) product between two vectors or matrices with identical sizes. The notation $\mathbf{I}_{N}$ is used to denote the $N\times N$ identity matrix. For set ${\mathcal{A}}$, $|{\mathcal{A}}|$ denotes its cardinality. For random variable $x$, ${\sf E}(x)$ and ${\sf Var}(x)$ denote the expectation and variance of $x$ respectively, while ${\sf E}(x\mid y)$ and ${\sf Var}(x\mid y)$ denote the expectation and variance of $x$ conditioned on $y$ respectively. The distribution of a proper complex Gaussian random variable $z$ with mean $\mu$ and variance $\nu$ is expressed by $z\sim\mathcal{CN}(z;\mu,\nu)=\frac{1}{\pi\nu}\exp{\left(|z-\mu|^{2}/{\nu}\right)}$. Lastly, the notation $:=$ is used to define a symbol as the expression on its right side.

II System Model

We consider the scenario with $K$ single-antenna users served by a large-scale array with $N\geq K$ antennas in the same time frequency resource. Specifically, we study the multi-user detector design. In the considered scenario, each user encodes its own information bit stream and modulates it to be a sequence of constellation points in $\mathcal{X}$ (e.g., 16-QAM). For a certain time slot, we denote the transmit symbol of the $k$-th user by $x_{k}$, and stack the transmit symbols of $K$ users by the vector $\mathbf{x}=[x_{1},x_{2},\ldots,x_{K}]^{T}$, where $\mathbf{x}\in\mathcal{X}^{K}$. The average energy of the transmit symbols of different users is assumed to be the same and denoted by $E_{x}$. The narrowband111Although the detector is proposed for narrowband channel, its use in wideband channel can be supported by utilizing orthogonal frequency division multiplexing (OFDM) technique, which decomposes the convolutive multipath channels into a bank of parallel flat-fading subchannels. equivalent baseband input-output relationship is given by

[TABLE]

where $\mathbf{y}\in\mathbb{C}^{N}$ denotes the signal vector received by the antenna array, $\mathbf{H}\in\mathbb{C}^{N\times K}$ is the channel matrix, and $\mathbf{n}\in\mathbb{C}^{N}$ represents the additive white Gaussian noise (AWGN) vector with distribution $\mathcal{CN}(\mathbf{0},\sigma^{2}\mathbf{I}_{N})$. We aim to produce a reliable estimate of $\mathbf{x}$ with the knowledge of received signals $\mathbf{y}$ and the channel matrix $\mathbf{H}$ for hard or soft decision.

A promising solution for dealing with the excessively large data load and computational complexity caused by the order-of-magnitude increase in the antenna array dimension is to partition the antenna array into a certain number of subarrays, each associated with an independent processing entity. Following this idea, we partition the whole array of $N$ antennas into $C\geq 1$ subarrays. Let $N_{c}$ denote the number of antennas in subarray $c\in\mathcal{C}=\{1,2,\ldots,C\}$, and thus $N=\sum_{c=1}^{C}N_{c}$. Correspondingly, the received and AWGN vectors and channel matrix can be partitioned as $\mathbf{y}=[\mathbf{y}_{1}^{T},\mathbf{y}_{2}^{T},\ldots,\mathbf{y}_{C}^{T}]^{T}$, $\mathbf{H}=[\mathbf{H}_{1}^{T},\mathbf{H}_{2}^{T},\ldots,\mathbf{H}_{C}^{T}]^{T}$, and $\mathbf{n}=[\mathbf{n}_{1}^{T},\mathbf{n}_{2}^{T},\ldots,\mathbf{n}_{C}^{T}]^{T}$. Then, the received vector corresponding to the $c$-th subarray is given by

[TABLE]

where ${\mathbf{H}_{c}\in\mathbb{C}^{N_{c}\times K}}$ denotes the local channel matrix between $K$ users and antennas of subarray $c$, and ${\mathbf{n}_{c}\in\mathbb{C}^{N_{c}}}$ is noise vector corresponding to the $c$-th subarray.

We call (2) the full subarray model in the subsequent sections because it involves full channel matrices. The sparsity in the local channel matrices inherent in the spatial non-stationarity can be exploited to further reduce complexity. When non-stationarity occurs, energy received from each user is concentrated on a small portion of subarrays, which is usually significantly smaller than $C$. Specifically, only a limited portion of users can be viewed as being served by a specific subarray. Therefore, the local channel matrices $\{\mathbf{H}_{c}\}_{c=1}^{C}$ are approximately column-wise sparse, that is, the columns of $\mathbf{H}_{c}$ corresponding the users with sufficiently large received energy on the $c$-th subarray have major values, while the other columns are close to zero. The local channel matrix $\mathbf{H}_{c}$ of each subarray can be estimated by applying compressive sensing techniques for multiple measurement vector problems (see [25, Section 11.6] and references therein) and then be trimmed by removing these columns with values close to zero. Accordingly, we denote the trimmed local channel matrix of subarray $c$ by $\tilde{\mathbf{H}}_{c}$, and we define $\mathcal{K}_{c}$ as the subset of index corresponding to the columns of $\mathbf{H}_{c}$ which are retained in $\tilde{\mathbf{H}}_{c}$. The matrix $\tilde{\mathbf{H}}_{c}$ is a $N_{c}\times K_{c}$ submatrix of $\mathbf{H}_{c}$ that comprises $K_{c}$ columns of $\mathbf{H}_{c}$ indexed by $\mathcal{K}_{c}$. We then let $K_{c}=|\mathcal{K}_{c}|$ and $\mathbf{x}_{c}=[\mathbf{x}]_{\mathcal{K}_{c}}$. Typically, $K_{c}$ is much smaller than $K$. With the above notations, the received model for each subarray is approximated by

[TABLE]

which is called as trimmed subarray model. The reduction in the dimensions of local channel matrices will lead to reduced complexity in signal detection.

The subarray-based baseband processing architecture is shown in Fig. 1. Each subarray is associated with its local processing module (LPM), which accesses only their local channel matrix and received vector. In particular, the $c$-th subarray produces a coarse estimate of $\mathbf{x}$ (or $\mathbf{x}_{c}$) based on $\mathbf{y}_{c}$ and $\mathbf{H}_{c}$ (or $\tilde{\mathbf{H}}_{c}$)222This paper mainly focuses on the detector design, thus we assume the local channel matrix, $\mathbf{H}_{c}$ or $\tilde{\mathbf{H}}_{c}$, is perfectly known at each subarray. In practical systems, the local channel matrix of each subarray can be acquired by uplink pilot [4]. which is feedforward to the central processing module (CPM). The CPM is responsible for generating consensus information based on the information delivered by each subarray (solid blue lines), which is then broadcast back to all subarrays (dash green lines), and producing the final output, which is then applied to generate the input of the channel decoder. The blue and green lines may correspond to dedicated data links, implemented by cables or optic fibers, or the data interfaces among chips or software modules depending on different architectures. In the subsequent sections, we focus on detector design based on this subarray-based architecture.

III EP Detector for Subarray-based Architecture

In this section, we exploit the EP principle to develop the detector for subarray-based architecture. We consider the fundamental full subarray model (2) to establish the technical foundation and facilitate a few theoretical analyses from the signal processing perspective. We start by providing detailed interpretations of the proposed detector, followed by intuitive convergence justification and fixed point characterization. Discussions related to efficient implementation of the proposed algorithm will be presented in the next section.

III-A Algorithm Derivation

We develop a detector by constructing a reliable estimate of $\mathbf{x}$ on the basis of the classical Bayesian inference framework [26], which is initiated by calculating the a posteriori distribution as follows:

[TABLE]

where $p(\mathbf{y}|\mathbf{x},\mathbf{H})$ denotes the likelihood function with the known channel matrix $\mathbf{H}$ and $p(\mathbf{x})$ represents the a priori distribution of $\mathbf{x}$. Without loss of generality, we consider the case where $\{x_{k}\}_{k=1}^{K}$ are drawn independently from $\mathcal{X}$ with equal probabilities, namely, $p(x_{k}=x)=1/|\mathcal{X}|$ for $\forall x\in\mathcal{X}$. We can then attain excellent inference in accordance with several optimal criteria. For example, we can achieve the minimum mean-square error (MMSE) criterion by computing the a posteriori distribution of $\mathbf{x}$. We can also achieve the maximum a posteriori (MAP) criterion by searching $\mathbf{x}$ over $\mathcal{X}^{N}$ which maximizes the a posteriori distribution in (4). However, direct computation of the aforementioned estimates is intractable because of the calculation of the high-dimensional integral or exponential complexity. This situation motivates us to pursue advanced mathematical tools to effectively approximate (4) with tractable complexity and appropriately fit the subarray-based architecture.

EP is an effective tool that can be used to approximate complicated distributions, such as (4), fitting appropriately the subarray-based architecture. The basic idea is to represent (4) using a graphical model called factor graph associated with a certain factorization of (4). Compared with the conventional EP-based MIMO detector [17], the factorization of (4) in this work further considers the subarray-based architecture (2), which is given by

[TABLE]

with ${f_{0}}(\mathbf{x})=p(\mathbf{x})$ and ${f_{c}}(\mathbf{x})=\exp{\left(-\|\mathbf{y}_{c}-\mathbf{H}_{c}\mathbf{x}\|^{2}/\sigma^{2}\right)}$ for $c\in\mathcal{C}$. In this manner, we represent (4) as a non-loopy factor graph with vector-valued nodes in Fig. 2a. This factor graph contains $C+1$ factor nodes $\{f_{0},f_{1},\ldots,f_{C}\}$ , and one variable node $\mathbf{x}$. Then, we restrict messages updated and transfered between nodes on the factor graph as Gaussian distributions. In principle, Gaussian distributions can be fully characterized by their expectations and variances. Thus, we only need to compute and propagate the expectations and variances of these messages. On this basis, this process is named as expectation propagation. Then, we approximate the a posteriori distribution (4) by computing and passing messages among the nodes in the factor graph (Fig. 2a) in an iterative manner based on certain rules, and we obtain Algorithm 1. We elaborate the derivation of Algorithm 1 in Appendix A, while provide in the next subsection some operational explanations for the proposed detector.

Remark 1

Various algorithms proposed previously can be derived utilizing the EP principle for certain forms of factor graph, in the same manner as the framework presented in Appendix A. For example, the procedure in Appendix A can be followed to derive the original centralized EP MIMO detector [17] with the factor graph shown in Fig. 2b by alternating between the selection of the factors $p(\mathbf{y}|\mathbf{x},\mathbf{H})$ and $p(\mathbf{x})$. In other words, when $C=1$, Algorithm 1 presented in the next subsection is reduced to the centralized EP MIMO detector [17]. Algorithms proposed in [22, 27, 23, 24] can be derived on the basis of the EP principle for a factor graph similar to Fig. 2b as well. By introducing an auxiliary variable and a factor node corresponding to the linear transformation, the EP principle can be exploited to develop effective algorithms involving nonlinear measurements, such as [28, 29, 30, 31]. Therefore, the EP principle will continue to help in developing effective algorithms for a great number of other applications, if we can describe the particular application scenario properly with a factor graph allowing for feasible computation of messages.

III-B Proposed Detector

The block diagram of the proposed EP-based detector is illustrated in Fig. 3, which provides detailed illustration for the block “Local Det.” in Fig. 1. The algorithm alternates between parallel computation in each LPM (Steps (1)$-$(3)), and combination and estimation in the CPM (Steps (4) and (5)). Each LPM in parallel produces coarse estimates of $\mathbf{x}$ based on their associated received vector, channel matrix, and feedback from the CPM of the last iteration. Then, the extrinsic information, which is passed forward to the CPM (shown by the blue solid arrow in Figs. 1 and 3), is calculated. Subsequently, the CPM refines the estimate of $\mathbf{x}$ considering the a priori distribution $p(\mathbf{x})$, which is broadcasted back to each LPM (shown by the green dotted arrow in Figs. 1 and 3) for the next iteration. The entire procedure is summarized mathematically in Algorithm 1 with explicit expressions for Steps (1)$-$(5). Then, we explain in detail the operations of the proposed algorithm.

We start by explaining the operations performed in the LPMs of each subarray. At the beginning of the iterative process, $\tau_{c}$ and $\boldsymbol{\gamma}_{c}$ are initialized to be $\mathbf{0}$ and $E^{-1}_{x}$ respectively. At this moment, no extra information about $\mathbf{x}$ is available except the mean and covariance of $\mathbf{x}$. Then, we compute the LMMSE estimation of $\mathbf{x}$ with respect to (w.r.t.) the linear model in (2) with the a priori mean and covariance matrix given by $\boldsymbol{\gamma}_{c}$ and $\tau^{-1}_{c}\mathbf{I}_{K}$ respectively. Following [26, Theorem 12.1] with some manipulations on the basis of the matrix inversion lemma, we can obtain the a posteriori covariance matrix and mean of $\mathbf{x}$ by (7a) and (7b), respectively. Next, we calculate the average a posteriori variance and take its reciprocal as $\omega_{c}$ in (7c). Subsequently, we calculate the extrinsic information, which can be obtained by excluding the a priori distribution of $\mathbf{x}$ from its a posteriori distribution. Following the Gaussian message combining rule [32, (54) and (55)], the extrinsic variance and mean can be given by (8a) and (8b) respectively. They are then transferred to the CPM via the feedforward link.

Subsequently, we turn to the operations performed by the CPM. After collecting all extrinsic means and variances from each subarray, the CPM performs the maximum-ratio combining (MRC) expressed by (9a) and (9b)333 The name “MRC” originates from the fact that Eqs. (9a) and (9b) are identical to the linear combination of \{\buildrel\lower 3.0pt\hbox{\scriptscriptstyle\rightarrow}\over{\mathbf{x}}_{c}\}_{c=1}^{C} with minimum post-combination signal-to-noise ratio (SNR) by viewing \buildrel\lower 3.0pt\hbox{\scriptscriptstyle\rightarrow}\over{\mathbf{x}}_{c} as AWGN observation of $\mathbf{x}$ with noise power $\eta^{-1}_{c}$.. The above MRC forms $\boldsymbol{\gamma}_{0}$ as an AWGN observation of $\mathbf{x}$ with noise power $\tau^{-1}_{0}$, that is

[TABLE]

where $\mathbf{n}_{0}\sim\mathcal{CN}(\mathbf{0},\tau^{-1}_{0}\mathbf{I})$. Together with the knowledge that each entry of $\mathbf{x}$ are equi-probably taken from the constellation $\mathcal{X}$, we prepare to calculate the MMSE estimation of $\mathbf{x}$ from (11). Equivalently, we compute the a posteriori mean and variance of $x_{k}$ for $k\in\mathcal{K}$ w.r.t. the posterior probability

[TABLE]

then the explicit expressions of $\hat{x}_{0,k}$ and $v_{0,k}$ can be given by

[TABLE]

respectively. Finally, the estimate $\hat{\mathbf{x}}_{0}$ and its corresponding average MSE computed in (10c) are sent back to the LPMs of each subarray through the feedback link. Next, each subarray calculates their a priori mean and variance of $\mathbf{x}$ as (6b) and (6a) respectively, thereby starting the next iteration. The above process is performed iteratively until a certain iteration stopping criterion is met. The most convenient stopping criterion is the maximum iteration numbers. After the convergence of the iteration, we obtain the approximated marginal posterior probability $\mathrm{P}(x_{k}\mid\mathbf{y},\mathbf{H})\approx\mathcal{CN}(x_{k};\hat{x}_{0,k},v_{0,k})$.

When computing \buildrel\lower 3.0pt\hbox{\scriptscriptstyle\rightarrow}\over{\mathbf{x}}_{c} in (8b), $\eta_{c}$ is divided. The subsequent operations in (9b) and (6b) include \eta_{c}\buildrel\lower 3.0pt\hbox{\scriptscriptstyle\rightarrow}\over{\mathbf{x}}_{c}. To avoid redundant division and multiplication of $\eta_{c}$, in practical applications, we can compute \buildrel\lower 3.0pt\hbox{\scriptscriptstyle\rightarrow}\over{\mathbf{x}}_{c} by \buildrel\lower 3.0pt\hbox{\scriptscriptstyle\rightarrow}\over{\mathbf{x}}_{c}=\omega_{c}\hat{\mathbf{x}}_{c}-\tau_{c}\boldsymbol{\gamma}_{c}. Accordingly, the computation of $\boldsymbol{\gamma}_{0}$ and $\boldsymbol{\gamma}_{c}$ in (8b) and (6b) can be simplified to be

[TABLE]

III-C Intuitive Verification of Convergence

Firstly, we characterize the performance for the proposed EP-based detector based on the evolution technique. Its basic idea is to select a few key parameters, which can be applied to characterize the statistical behavior of the iterative process and then to calculate the so-called transfer functions to track their evolution via iterations. We can observe from Algorithm 1 that the final estimate of $\mathbf{x}$ is produced from their AWGN observation $\boldsymbol{\gamma}_{0}$ and that its performance can be determined by the average noise power $\tau^{-1}_{0}$ from well-established formulas. Therefore, the first parameter is selected as $\rho=\tau^{-1}_{0}$. In addition, from (8a) and (9a), we find that $\tau_{0}$ and $\{\tau_{c}\}_{c=1}^{C}$ are mutually related in a recursive manner. Hence they should also be selected as the parameters to be examined. To enable concise expression, we denote $\nu_{c}=\tau^{-1}_{c}$ for $c\in\mathcal{C}$. Subsequently, we calculate transfer functions $\rho$ and $\nu_{c}$’s to show their evolution with the iteration, thereby inducing the following proposition.

Lemma 1

The evolution of parameters $\rho$ and $\nu_{c}$’s can be characterized recursively by

[TABLE]

which is initiated by $\nu^{1}_{c}=E_{x}$ for $\forall c\in\mathcal{C}$, where the subscript $t$ denotes the iteration index, $\lambda_{i,c}$ is the $i$th eigenvalue of the matrix $\mathbf{H}^{H}_{c}\mathbf{H}_{c}$, ${\sf mse}_{c}\left(\nu^{t}_{c}\right)$ denotes the average MSE of the estimator (7c) given $\mathbf{H}_{c}$ and $\nu^{t}_{c}$ for each subarray $c\in\mathcal{C}$, which is given by

[TABLE]

and ${\sf mse}_{0}\left(\rho^{t}\right)$ denotes the average MSE of the MMSE estimation of $x$ given its AWGN observation $r$ with average noise power $\rho^{t}$, which is given by

[TABLE]

Proof:

The transfer function of $\rho$ from $\nu_{c}$’s in (13a) can be calculated straightforwardly by substituting (8a) into (9a). Note that here we compute the average a posteriori variance $\omega_{c}$ of the estimator (7c) by applying the singular-value decomposition (SVD) of $\mathbf{H}_{c}$. The transfer functions of $\nu_{c}$’s from $\rho$ can be calculated similarly by substituting (8a) into (6a). For a sufficiently large $N$, the average a posteriori variance $\omega_{0}$ in (10c) can be computed compactly by (15) from the central limit theorem. ∎

The above evolution analysis has many interesting applications, such as intuitive convergence verification. Proving the overall convergence of the proposed EP-based detector is generally difficult [19]. However, the result of the evolution analysis provides us with the possibility to justify the convergence intuitively if we can show that the sequences $\{\rho^{t}\}$ and $\{\nu^{t}_{c}\}_{c=1}^{C}$ w.r.t. $t$ are bounded and monotonous, as shown in the proposition below. Later in Section V, we will show through simulation results that the evolution analysis can provide precise performance prediction, and that the proposed algorithm converges.

Proposition 1

The sequences $\{\rho^{t}\}$ and $\{\nu^{t}_{c}\}_{c=1}^{C}$ w.r.t. $t$ are bounded and monotonically decreasing.

Proof:

See Appendix B. ∎

III-D Fixed Point Characterization

In this subsection, we characterize the fixed points of the proposed EP-based detector. In particular, we show that the fixed points are stationary points of a relaxed version of the KL divergence minimization distribution approximation problem. We then clarify the relation of the fixed points to the replica prediction of the asymptotic MMSE performance. Before proceeding, a simple consistency result, which is needed by subsequent discussions, is given by Lemma 2.

Lemma 2

For any fixed point of Algorithm 1 with $\tau_{0}+\sum_{c=1}^{C}\tau_{c}>0$, we have

[TABLE]

Proof:

Substituting (8a) into (6a), we have $\omega_{0}=\omega_{c}$ for $c\in\mathcal{C}$. Similarly, substituting (8b) into (6b), and together with the fact that $\omega_{0}=\omega_{c}$, we have $\hat{\mathbf{x}}_{0}=\hat{\mathbf{x}}_{c}$ for $c\in\mathcal{C}$. Then, (16a) and (16b) can be proved by combining (8a) with (9a), and combining (8b) with (9b) respectively. ∎

The typical task of approximating a complicated distribution, such as (4), is to minimize the Kullback-Leibler (KL) divergence between the approximated distribution and the original one over a given distribution family. However, this minimization is generally intractable as it involves a search over a family of $K$-dimension distributions. However, exploiting certain relaxations of the aforementioned KL divergence minimization problem allows for feasible solutions. The following proposition reveals that the fixed points of the proposed detector characterizes the stationary point of a particular relaxed version of the KL divergence minimization problem.

Proposition 2

Denote the probability density functions $b_{0}(\mathbf{x})$, $\{b_{c}(\mathbf{x})\}_{c=1}^{C}$, and $q(\mathbf{x})$ parameterized by the common values $\omega$ and $\hat{\mathbf{x}}$ of the fixed points shown in (16a) and (16b) as

[TABLE]

respectively, where $Z_{0}\left(\boldsymbol{\gamma}_{0}\right)$, $\{Z_{c}\left(\boldsymbol{\gamma}_{c}\right)\}_{c=1}^{C}$, and $Z_{q}\left(\hat{\mathbf{x}}\right)$ are the normalization factor for their corresponding density functions. Then, $b_{0}(\mathbf{x})$, $\{b_{c}(\mathbf{x})\}_{c=1}^{C}$, and $q(\mathbf{x})$ are the stationary points of the optimization of the Bethe free energy subject to the moment-matching constraints, which can be expressed as

[TABLE]

where equality constraints (18b) and (18c) are set for $c=0,1,\ldots,C$, ${\sf KL}\left({\cdot\left\|\cdot\right.}\right)$ denotes the KL divergence of two distributions and $H(\cdot)$ denotes the differential entropy. In addition, $\omega$ and $\hat{\mathbf{x}}$ satisfy that

[TABLE]

Proof:

See Appendix C. ∎

Now, we show how the KL divergence minimization problem is relaxed to be problem (18). Under the constraint that $b_{1}(\mathbf{x})=\cdots=b_{C}(\mathbf{x})=q(\mathbf{x})$, the minimization of $J\left(b_{1},\ldots,b_{C},q\right)$ in (18a) is equivalent to approximate $p(\mathbf{x}|\mathbf{y},\mathbf{H})$ in (4) by a distribution with minimum KL divergence with $p(\mathbf{x}|\mathbf{y},\mathbf{H})$, which is still intractable. We further relax the constraint $b_{1}(\mathbf{x})=\cdots=b_{C}(\mathbf{x})=q(\mathbf{x})$ to be the moment-matching constraint in (18b) and (18c), which requires a match in their first moments and average in their second moments. Then, we provide an intuitive explanation of (11) as follows.

Remark 2

Consider the belief estimate $b_{0}(\mathbf{x})$ given by (17a). If $\boldsymbol{\gamma}_{0}$ is modeled as a random vector, then (17a) implies that $b_{0}(\mathbf{x})$ represents the a posteriori of $\mathbf{x}$ given $\boldsymbol{\gamma}_{0}$, namely, $p\left(\mathbf{x}|\boldsymbol{\gamma}_{0}\right)$. Then, from the Bayes rule, we can express the likelihood function as

[TABLE]

where $p\left(\boldsymbol{\gamma}_{0}\right)$ denotes the marginal distribution $p\left(\boldsymbol{\gamma}_{0}\right)=\int p\left(\mathbf{x}|\boldsymbol{\gamma}_{0}\right)p\left(\mathbf{x}\right)\mathrm{d}\mathbf{x}$. An admissible choice of $p\left(\boldsymbol{\gamma}_{0}\right)$ that satisfies $p\left(\boldsymbol{\gamma}_{0}\right)\geq 0$, $\int p\left(\boldsymbol{\gamma}_{0}\right)\mathrm{d}\boldsymbol{\gamma}_{0}=1$, and $\int p\left(\boldsymbol{\gamma}_{0}|\mathbf{x}\right)\mathrm{d}\boldsymbol{\gamma}_{0}=1$ is $p\left(\boldsymbol{\gamma}_{0}\right)=\pi^{-N}\tau^{N}_{0}Z_{0}\left(\boldsymbol{\gamma}_{0}\right)$. Then, we have $p\left(\boldsymbol{\gamma}_{0}|\mathbf{x}\right)=\mathcal{CN}\left(\boldsymbol{\gamma}_{0};\mathbf{x},\tau_{0}\mathbf{I}\right)$. In this case, we can interpret $\boldsymbol{\gamma}_{0}$ as an unbiased estimate of $\mathbf{x}$ Gaussian estimation error of variance $\tau_{0}$. On this basis, we can calculate the log-likelihood ratio (LLR) through $\boldsymbol{\gamma}_{0}$ and $\tau_{0}$ for soft decoding. In addition, the hard decision can also be performed by determining $x_{k}\in\mathcal{X}$ for each $k\in\mathcal{K}$ with the shortest distance to $\gamma_{0,k}$, as we do in Algorithm 1.

Next, we characterize the asymptotic performance of the proposed EP-based detector as the following proposition.

Proposition 3

When $\mathbf{H}$ is a random matrix with i.i.d. elements, for any fixed point $\omega$ and $\tau_{0}$ of proposed EP-based detector, it holds that

[TABLE]

under the large system regime where $N,N_{1},N_{2},\ldots,N_{C},K\to\infty$ with $K/N$ being a fixed constant.

Proof:

See Appendix D. ∎

Remark 3

*We find that (20) is identical to the fixed-point equation characterization of the asymptotic MSE performance for (4) derived from the replica method [33, Eq. (17)]. Therefore, when the channel matrix $\mathbf{H}$ contains i.i.d. entries, the proposed EP detector can potentially achieve the MMSE performance in certain asymptotic and random regimes. In other words, under any number and size of subarray, the proposed detector will produce similar performance. Furthermore, simulation results show that even when $\mathbf{H}$ is non-i.i.d., reducing $N_{c}$ does not result in significant performance loss. As discussed in Section IV-B, including a small number of antennas in each subarray allows for simple-form computation of some matrix operations. The above observations inspire us to parallelize the complex operation into several simple-form computations, giving rise to the hierarchical implementation architecture proposed in Section IV-B. *

IV Efficient implementation-related issues

In this section, we shift our discussion to the engineering perspective. In Subsection IV-A, we demonstrate that the computational complexity and the amount of data required to be transfered from each subarrays to the CPM can be reduced by using the trimmed local channel matrices. Then, in Subsection IV-B, we propose a hierarchical architecture that enables the parallelism improvement. In Subsection IV-C, we modify the proposed EP-based detector so that only one feedforward from subarrays to the CPM is required in the entire iteration procedure. Finally in Subsection IV-D, we provide complexity comparison.

IV-A Implementation of Algorithm 1 based on trimmed subarray model

The principal advantage of implementing Algorithm 1 based on the trimmed subarray model (3) is the reduction in computational complexity and the amount of data needed to transfer from each subarray to the CPM. In this case, Step (2) in Algorithm 1 produces $\hat{\mathbf{x}}_{c}$ as the LMMSE estimate of $\mathbf{x}_{c}$, which is given by

[TABLE]

Then, the number of complex-value multiplication required in producing $\hat{\mathbf{x}}_{c}$ is reduced from $K(N_{c}+K)$ to $K_{c}(N_{c}+K_{c})$, and the number of complex-value multiplication required in computing $\boldsymbol{\Sigma}_{c}$ using the method presented in Algorithm 2 is reduced from $N_{c}K(2K+1)$ to $N_{c}K_{c}(2K_{c}+1)$. Moreover, the number of complex numbers required to transfer to the CPM, namely the dimension of \buildrel\lower 3.0pt\hbox{\scriptscriptstyle\rightarrow}\over{\mathbf{x}}_{c}, is reduced from $K$ to $K_{c}$. Typically, $K_{c}$ is much smaller than $K$. Therefore, the computational complexity and information exchange overhead can be considerably reduced. Accordingly, the MRC combination performed in Step (4) of Algorithm 1 needs to be modified as

[TABLE]

for $k\in\mathcal{K}$, where $\mathcal{C}_{k}$ denotes the subset of $\mathcal{C}$ satisfying the condition that if $c\in\mathcal{C}_{k}$, then $k\in\mathcal{K}_{c}$. In the computations of $\hat{x}_{0,k}$ and $v_{0,k}$ using (12), we substitute $\tau_{0}$ by $\tau_{0,k}$. In addition, we compute $\boldsymbol{\gamma}_{c}$ as \boldsymbol{\gamma}_{c}=\tau^{-1}_{c}\left(\omega_{0}[\hat{\mathbf{x}}_{0}]_{\mathcal{K}_{c}}-\eta_{c}\buildrel\lower 3.0pt\hbox{\scriptscriptstyle\rightarrow}\over{\mathbf{x}}_{c}\right). Except for the steps mentioned above, the expressions for the other steps are unchanged.

IV-B Hierarchical Implementation Architecture

The primary bottleneck in the practical implementation of the proposed EP-based detector is that the matrix inversion in (7a), whose complexity is in the order of $K^{3}$ considering full subarray model, or $K_{c}^{3}$ considering trimmed subarray model, is required in every iteration. Under a crowded scenario, this complexity is unaffordable. However, when the antenna number in each subarray is small, a simple computation of the matrix inversion can be performed. Consider the extreme case, $N_{c}=1$, in which a computationally efficient formula is available for the matrix inversion. In this case, $\mathbf{H}_{c}$ is the $c$-th row of $\mathbf{H}$. The matrix $\mathbf{H}^{H}$ can be blocked as $\mathbf{H}^{H}=[\bar{\mathbf{h}}_{1},\bar{\mathbf{h}}_{2},\ldots,\bar{\mathbf{h}}_{C}]$. Then $\boldsymbol{\Sigma}_{c}$ in (7a) can be calculated in a very simple form shown below,

[TABLE]

which involves the multiplication of a column and a row vector divided by a scalar. When $N_{c}$ is adequately small444What number of $N_{c}$ can be called adequately small depends on the computational speed of the component in which the algorithm is implemented., the matrix inversion can be computed by recursively using the above formula for $N_{c}$ times by exploiting the following observation:

[TABLE]

We let $\mathbf{A}_{c,0}=\tau_{c}\mathbf{I}_{K}$ and $\mathbf{A}_{c,j}=\tau_{c}\mathbf{I}_{K}+\sigma^{-2}\sum_{i=1}^{j}\bar{\mathbf{h}}_{i,c}\bar{\mathbf{h}}^{H}_{i,c}$. Evidently, $\boldsymbol{\Sigma}_{c}=\mathbf{A}^{-1}_{c,N_{c}}$ and $\mathbf{A}_{c,j}=\mathbf{A}_{c,j-1}+\sigma^{-2}\bar{\mathbf{h}}_{j,c}\bar{\mathbf{h}}^{H}_{j,c}$. If $\mathbf{A}^{-1}_{c,j-1}$ is given, $\mathbf{A}^{-1}_{c,j}$ can be then obtained in the same manner as (22). On the basis of the above observation, $\boldsymbol{\Sigma}_{c}$ can be calculated recursively using the procedure given by Algorithm 2.

If $N_{c}$ is a very small number, such as 1, 2, or 4, and $N$ is relatively large, such as 512, the number of subarrays is excessive. In this condition, the amount of data that needs to be interchanged through the interface between subarrays and the CPM will be large. In some cases, $N_{c}$ is fixed by hardware constraint, which cannot be adapted flexibly. Then, the above strategy to make $N_{c}$ small is not feasible. To keep $C$ a reasonable number but let the above strategy work, we propose to implement the proposed algorithm in a hierarchical manner by further dividing each subarray into several secondary subarrays. For example, we divide $N_{c}$ antennas in the $c$-th subarray into $M_{c}$ secondary subarrays. Then, we denote the antenna number in the $c^{\prime}$-th secondary subarray by $N_{c,c^{\prime}}$. Although $N_{c}$ may be a relatively large number, we can set $N_{c,c^{\prime}}$ to be small numbers such as 1, 2, or 4. The LPM of each subarray assigns $M_{c}$ parallel computation units for each secondary subarray. Then $\mathbf{y}_{c}$ and $\mathbf{H}_{c}$ can be further blocked as $\mathbf{y}_{c}=[\mathbf{y}_{c,1}^{T},\mathbf{y}_{c,2}^{T},\ldots,\mathbf{y}_{c,M_{c}}^{T}]^{T}$ and $\mathbf{H}_{c}=[\mathbf{H}_{c,1}^{T},\mathbf{H}_{c,2}^{T},\ldots,\mathbf{H}_{c,M_{c}}^{T}]^{T}$ respectively. In each iteration, each secondary subarray performs Steps (1) - (3) of Algorithm 1, substituting the subscripts of all variables with subscript ’c’ by ’${}_{c,c^{\prime}}$’. After the extrinsic information of each secondary subarray $\eta_{c,c^{\prime}}$ and \buildrel\lower 3.0pt\hbox{\scriptscriptstyle\rightarrow}\over{\mathbf{x}}_{c,c^{\prime}} are obtained, each subarray combines them as its extrinsic information to deliver to the CPM as

[TABLE]

After the feedback of $\omega_{0}$ and $\hat{\mathbf{x}}_{0}$ from the CPM, each subarray broadcasts them to its secondary subarrays for the next iteration. The above hierarchical operation is equivalent to keeping a small number $N_{c}$ of antennas in each subarray. Notably, the above discussion is for the algorithm corresponding to the full subarray model and the same way can be directly applied for the algorithm for trimmed subarray model.

IV-C Modified EP-based detector requiring one feedforward

Another practical challenge is the need of multiple feedforward ($\{\eta_{c}\}_{c=1}^{C}$ and \{\buildrel\lower 3.0pt\hbox{\scriptscriptstyle\rightarrow}\over{\mathbf{x}}_{c}\}_{c=1}^{C} from each subarray to the CPM) and feedback ($\omega_{0}$ and $\hat{\mathbf{x}}_{0}$ from the CPM to each subarray) in every iteration for the information exchange between the subarrays and the CPM. However, for application scenarios where the hardware entities of the subarrays and the CPM are placed distantly and connected by dedicated fronthaul links, such frequent iterative information exchange suffers from high interconnect latency. To address this challenge, we modify the proposed EP-based detector so that the entire iterative procedure requires only one feedforward. To achieve this aim, two potential schemes are presented as follows.

• •

For each received signal $\mathbf{y}$, we perform only one iteration of Algorithm 1 to obtain the final output.

• •

Each subarray performs the entire procedure of Algorithm 1 (or the centralized EP algorithm [17]) for their local input-output relationship $\mathbf{y}_{c}=\tilde{\mathbf{H}}_{c}\mathbf{x}_{c}+\mathbf{n}_{c}$ in parallel. After convergence, each subarray transfers the estimate of $\mathbf{x}$ and the corresponding variance (corresponding to the “output” line in Algorithm 1) to the CPM for MRC combination (in the same form of Eqs. (9a) and (9b)) as the final output.

The first scheme is the most straightforward idea to achieve one feedforward. However, it prevents performance improvement via iterations inherently provided by EP-based algorithms, which results in performance degradation. This observation motivates us to retain the iteration and propose the second scheme. The drawback of this scheme is that the information exchange among subarrays is disabled. Notably, when non-stationarity occurs, each subarray serves only a small portion of active users, and the overlap among sets of users served by different subarrays is considerably small. Therefore, disabling information among subarrays may not result in a significant performance loss.

IV-D Implementation Complexity

Table I provides quantitative results in assessing implementation complexity of the proposed subarray-based detector. We compare the complexity of Algorithm 1, one feedforward architecture proposed in Section IV-C, and decentralized EP detector proposed in [17]. The complexity of the operations performed in LPM and CPM are counted separately. Operations that implement in a single LPM are counted because all LPMs simultaneously perform parallel computations. Table I lists the total number of real-valued multiplications, exponential operations, and real numbers exchanging between LPMs and the CPM (denoted by Mult., Exp. and Trans., respectively) required for performing $T$ iterations. We also compare the complexity when full channel matrix $\{\mathbf{H}_{c}\}_{c=1}^{C}$ and trimmed channel matrix $\{\tilde{\mathbf{H}}_{c}\}_{c=1}^{C}$ are utilized by each LPM. Table I shows that the required number of multiplications in each LPM is proportional to $K^{2}$ or $K^{2}_{c}$, and the dimension of data required to transfer from each LPM to the CPM is proportional to $K$ or $K_{c}$. Due to the non-stationarity, $K_{c}$ is typically considerably smaller than $K$. Hence, the sparsity in local channel matrices reduces computational complexity and information exchange overhead. In addition, one feedforward architecture significantly reduces the information exchange overhead compared with Algorithm 1, thereby reducing the latency caused by frequent information exchange. Therefore, the one feedforward architecture is suitable for scenarios where the hardware entities corresponding to different subarrays and the CPM are distributed over a large area and interconnected by capacity limited fronthaul links.

V Simulation Results

In this section, we demonstrate numerical results for performance evaluation and comparison. We take the uncoded bit error rate (BER) as the performance metric using Monte Carlo simulations of 10,000 independent channel realizations. Without the loss of generality, we set the number of antennas in each subarray $N_{c}$ to be identical in our simulations. In addition, we mark the iteration numbers in the figures below to provide clues about the computational delays of each algorithm.

V-A Stationary Case

In this subsection, we simulate the basic stationary case where a moderate number of antennas are compactly installed. Various algorithms proposed for the full subarray model (2) are evaluated and compared to examine their lossless performance and verify the validity of our analysis. The stationary case corresponds to compact antenna deployment. Hence, the correlation among antennas are considered. Since stationary case corresponds to compact antenna deployment, we will consider correlation among antennas here. We generate the channel matrix in our simulations as $\mathbf{H}=\boldsymbol{\Sigma}_{R}^{\frac{1}{2}}\mathbf{H}_{R}$, where $\mathbf{H}_{R}$ is the Rayleigh fading matrix with i.i.d. Gaussian entries with zero mean and variance $K^{-1}$, and $\boldsymbol{\Sigma}_{R}$ characterizes the correlation among the elements of the antenna array, whose $(i,j)$-th entry is set to be $\kappa^{|i-j|}$, where $\kappa$ will be called correlation coefficients subsequently. Moreover, we set the number of antennas $N=64$ and the number of users $K=16$ and elements of $\mathbf{x}$ are drawn independently from the equiprobable 16-QAM constellation.

We first present the uncoded BER provided by Algorithm 1 versus the iteration numbers under different numbers of $N_{c}$ in Figs. 4a - 4d. The results of Rayleigh channel ($\kappa=0$) and correlated channel ($\kappa=0.5$) are shown. We observe from these figures that in all cases, the values of BER initially reduce and then remain stable with only minor fluctuations. This observation alleviates our worries about the possible numerical instability of the proposed algorithm with a small number of antennas in each subarray ($N_{c}=$1,2 and 4) and the correlated channel. The results in Fig. 4 justify our analysis in Section III-B. From Fig. 4, we also find that, with smaller number of antennas in each subarray, one or two more iterations are required to attain the numerical stability. Specifically, when $N_{c}=$ 1, 2, or 4, the numerical stability is attained with 5 or 6 iterations, whereas with larger $N_{c}=$ 16, 4 iterations are enough to attain the numerical stability. The decrease in BER is minor in the last few iterations. Therefore, we recommend performing 3 or 4 iterations rather than waiting until stability is achieved completely in the practical applications.

To examine the impacts of the value of $N_{c}$ on the performance of one iteration and that after stability, we provide BER versus $N_{c}$ with different iteration numbers in Fig. 5, where we include the performance of seven iterations as that after stability. As implied by Fig. 5a, larger $N_{c}$ leads to better BER performance when one iteration is performed. It finding is natural because one iteration of Algorithm 1 is equivalent to conducting MRC for the LMMSE detection for each subarray without any information exchange among different clusters, thereby improving the performance. Then, in the application scenario where only one iteration is allowed, having more antennas in each subarray is favorable. However, different results can be observed for the performance after stability. For the Rayleigh channel, the partition into subarrays does not result in performance loss compared with centralized case, namely, $N_{c}=64$, thereby justifying the results revealed in Proposition 3. Meanwhile, in the correlated channel, small $N_{c}$ values give rise to slight performance loss compared with larger value of $N_{c}$. We can conclude from the above observations that the computational complexity of matrix inversion can be alleviated by partitioning into subarrays with a small number of antennas in each of them at a low cost of performance degradation.

Then, in Fig. 6, we compare the performance of Algorithm 1 with that of various algorithms under full subarray model, including decentralized conjugate gradient (DeCG) and alternating direction method of multipliers (DeADMM) proposed in [34]555In our simulation, we used part of the code downloaded from https://github.com/VIP-Group/DBP, and the original centralized EP MIMO detector proposed in [17]. To facilitate parallel simple computation of matrix inversion using Algorithm 2, we set $N_{c}=2$, which yields a similar performance as that under larger $N_{c}$ as shown in Fig. 5b. In particular, the performances achieved with one iteration are emphasized because performing one iteration takes up the least computational resource and requires only one feedforward. We find from Fig. 6 that with one iteration, Algorithm 1 performs better than DeGC and DeADMM. Similarly, after attaining numerical stability, Algorithm 1 also outperforms DeGC and DeADMM. Moreover, we observe that the performance of Algorithm 1 approaches that of the centralized EP MIMO detector. This observation is consistent with the results shown in Fig. 5. Finally, when implementing DeADMM, we tend to choose larger $N_{c}$, as we observe that when $N_{c}=2$, DeADMM requires several times number of iterations to attain the same performance as that when $N_{c}=16$. Therefore, choosing smaller $N_{c}$ gives rise to further computational delay.

V-B Non-stationary Case

In this subsection, we consider a linear antenna array for our simulation setup, where the length of array is 250 meters with $N=512$ antennas, and $K=16$ active users have the same vertical distance to the array666Here, we follow the simulation setup in [12] to guarantee fair performance comparison with the algorithm proposed in [12]. and are uniformly distributed along the array. Under this setting, the received energy from each user varies among different portions of the antenna array, and spatial non-stationarity occurs. We then compare the performances of the algorithms that are proposed for the trimmed subarray model. On this basis, we generate the channel matrix as $\mathbf{H}=\mathbf{D}\odot\mathbf{H}_{R}$, where $\mathbf{H}_{R}$ is the Rayleigh fading matrix, and $\mathbf{D}$ contains the large-scale fading factors whose $(i,j)$-th entry is given by $d^{-1}_{ij}$ with $d_{ij}$ denoting the distance between $i$-th antenna and $j$-th user. The process of channel matrix acquisition is not simulated, therefore, we utilize [12, Algorithm 2] for the full channel matrix $\mathbf{H}$ with the power threshold set to 0.9 to determine which part of users is served by each subarray. In this manner, we generate trimmed channel matrices $\{\tilde{\mathbf{H}}_{c}\}$ used in our simulations.

Fig. 7 shows the BER versus the iteration numbers for the EP-based detector derived for the trimmed subarray model (Section IV-A) under different $N_{c}$ and SNRs. In all cases, the reductions in the BER values are minor after two iterations, meaning that numerical stability is reached. Thus, in the practical application, we recommend to perform two or three iterations. In addition, as the value of $N_{c}$ decreases, slight performance loss can be observed (this finding can also be observed in Fig. 8a). Hence, in practical scenario with non-i.i.d. channel matrix, implementing the hierarchical implementation architecture proposed in Section IV-B is only at a minor cost of BER performance.

As shown in Fig. 8, we compare the BER performances of different algorithms that are proposed for the trimmed subarray model, including the EP-based detector derived for the trimmed subarray model (Section IV-A), the EP-based and approximate message passing (AMP)-based one feedforward architecture 777These schemes refer to implement EP and AMP in parallel in each LPM by AMP in the one feedforward architecture presented in Section IV-C., and the SIC-based method proposed in [12, Algorithm 3]. Fig. 8a shows the curve of BER versus $N_{c}$, and Fig. 8b shows that of BER versus SNR. In general, the EP-based algorithms outperform AMP- and SIC-based algorithms proposed in [12, Algorithm 3]. When applying EP-based algorithms, we prefer to set a smaller $N_{c}$ because doing so allows us to parallelize the matrix inversion to reduce complexity. However, when applying AMP-based one feedforward method [12, Algorithm 3], we tend to set larger $N_{c}$ because setting smaller $N_{c}$ suffers from significant performance loss. In addition, in the small SNR regime, the performance degradation of the modified EP-based detector requiring one feedforward is minor. As SNR increases, this performance degradation further increases. Therefore, the one feedforward method is favorable in the small SNR regime.

VI Conclusion

In this paper, the EP principle was exploited for the derivation of an efficient detector of extra-large-scale massive MIMO systems with the subarray-based processing architecture. We described the a posteriori distribution as a factor graph and developed the iterative algorithm by computing and transferring messages among different nodes on the factor graph. We intuitively verified the convergence of the proposed algorithm via the evolution analysis. We also characterized the fixed points of the algorithm. We proved that these fixed points are the stationary points of the Bethe free energy optimization subject to the moment-matching constraints and are identical to the fixed point equations derived from the replica method. We also discussed many implementation-related aspects, such as how to exploit non-stationarity for complexity reduction, the hierarchical implementation architecture to allow for parallel simple computation of matrix inversion, and the design of an algorithm with only one feedforward. Finally, we demonstrated through simulation results that the proposed detector outperforms its counterparts and verified the validity of our analysis.

Appendix A Message-Passing Derivation of Algorithm 1

In this appendix, we present the message-passing derivation of Algorithm 1. Following the EP principle [20, 19, 21], we constrain the messages to dwell in the family ${\mathcal{F}}$ of Gaussian-distributed vectors with independent elements, which is in the form of $\mathcal{CN}\left(\mathbf{x};\mathbf{m},{\sf diag}(\mathbf{v})\right)$. For a factor graph with factor nodes $\{f_{0},f_{1},\ldots,f_{C}\}$, all connected to a set of variable node $\{\mathbf{x}_{i}\}$, the messages are computed and updated iteratively by performing the steps below.

Factor selection: One or a set of factors is selected from $\{f_{0},f_{1},\ldots,f_{C}\}$. 2. 2.

Variable-to-factor messages: The message(s) from a variable node (e.g., $\mathbf{x}_{j}$) to the selected factor node(s) (e.g., $f_{i}$) is computed by

[TABLE] 3. 3.

Factor-to-variable messages: To compute the message(s) from the selected factor node(s) (e.g., $f_{i}$) to a variable node (e.g., $\mathbf{x}_{j}$), the approximate belief $b_{f_{i}}(\mathbf{x}_{j})$ is initially computed at $f_{i}$, which is given by

[TABLE]

where the divergence measure for two distributions, such as $p(x)$ and $q(x)$, ${D_{1}}\left[{p\left(x\right)\left\|{q\left(x\right)}\right.}\right]$ is defined by

[TABLE]

Then, we set the message ${\mu_{{f_{i}}\to\mathbf{x}_{j}}}\left(\mathbf{x}_{j}\right)$ as ${\mu_{{f_{i}}\to\mathbf{x}_{j}}}\left(\mathbf{x}_{j}\right)\propto b_{f_{i}}(\mathbf{x}_{j})/{\mu_{\mathbf{x}_{j}\to{f_{i}}}(\mathbf{x}_{j})}$.

Remark 4

The factor graph in Fig. 2a has only one variable node $\mathbf{x}$, then (24) reduces to

[TABLE]

Now, we consider the minimization in (25). Following the result in [19], the belief $b_{f_{i}}(\mathbf{x})\in{\mathcal{F}}$ minimizing $D_{1}$ in (25) is obtained as $\mathcal{CN}\left(\mathbf{x};\mathbf{m}_{q},{\sf diag}(\mathbf{v}_{q})\right)$, where $\mathbf{m}_{q}$ and $\mathbf{v}_{q}$ denote the vectors comprising the mean and variance of each element of $\mathbf{x}$ w.r.t. the distribution ${{\mu_{\mathbf{x}\to{f_{i}}}}\left(\mathbf{x}\right){f_{i}}\left(\mathbf{x}\right)}$ respectively. In addition, similar to strategies utilized in [22, 27, 23, 24], we take the average of $\mathbf{v}_{q}$ as the variance for the belief, which yield $b_{f_{i}}(\mathbf{x})\propto\mathcal{CN}\left(\mathbf{x};\mathbf{m}_{q},K^{-1}{\sf tr}({\sf diag}(\mathbf{v}_{q}))\mathbf{I}_{K}\right)$, to save the computational resources when implementing the algorithms.

Then, we elaborate how to construct Algorithm 1 from the aforementioned rules. In accordance to the decentralized processing architecture in Fig. 1, we alternate the associated Steps (2) and (3) between the selection of factors $\{f_{1},\ldots,f_{C}\}$ and that of factor $f_{0}$, which corresponds to the processing in each LPM and that at the CPM. At the first iteration, we skip the Step (2) and initialize the message from the variable node $\mathbf{x}$ to factors $\{f_{1},f_{2},\ldots,f_{C}\}$ by ${\mu_{\mathbf{x}\to{f_{c}}}}\left(\mathbf{x}\right)=\mathcal{CN}\left(\mathbf{x};\mathbf{0},E_{x}\mathbf{I}\right)$ for $c\in\mathcal{C}$. Notably, in Algorithm 1, $\tau_{c}$ and $\boldsymbol{\gamma}_{c}$ denote the reciprocal of the variance and mean of messages $\{{\mu_{\mathbf{x}\to{f_{c}}}}\}_{c=0}^{C}$ respectively, namely, ${\mu_{\mathbf{x}\to{f_{c}}}}(\mathbf{x})=\mathcal{CN}\left(\mathbf{x};\boldsymbol{\gamma}_{c},\tau^{-1}_{c}\mathbf{I}\right)$. Therefore, the initialization of Algorithm 1 along with (6a) and (6b) yield the above initializations for ${\mu_{\mathbf{x}\to{f_{c}}}}\left(\mathbf{x}\right)$’s. Then, we proceed to Step (3), which starts by computing the approximate belief $\{b_{f_{c}}(\mathbf{x})\}_{c=0}^{C}$ for $c\in\mathcal{C}$. To solve the $D_{1}$-minimization problem in (23), we compute the a posteriori mean and covariance matrix w.r.t. the likelihood functions ${f_{c}}\left(\mathbf{x}\right)=\mathcal{CN}(\mathbf{y};\mathbf{H}_{c}\mathbf{x},\sigma^{2}\mathbf{I})$ and the Gaussian a priori ${\mu_{\mathbf{x}\to{f_{c}}}}\left(\mathbf{x}\right)$. From [26, Theorem 12.1], this a posteriori distribution is Gaussian with covariance matrix $\boldsymbol{\Sigma}_{c}$ and mean $\hat{\mathbf{x}}_{c}$ derived as (7a) and (7b), respectively. Then we take the reciprocal of the average diagonal elements for $\boldsymbol{\Sigma}_{c}$ as $\omega_{c}$ in (7c). Finally, we have $b_{f_{c}}(\mathbf{x})=\mathcal{CN}\left(\mathbf{x};\hat{\mathbf{x}}_{c},\omega^{-1}_{c}\mathbf{I}\right)$. Next, we set the message ${\mu_{{f_{c}}\to\mathbf{x}}}\left(\mathbf{x}\right)$ to be the ratio of the approximate belief $b_{f_{c}}(\mathbf{x})$ to the most recent message ${\mu_{\mathbf{x}\to{f_{c}}}}\left(\mathbf{x}\right)$, which is given by

[TABLE]

Then, we have {\mu_{{f_{c}}\to\mathbf{x}}}\left(\mathbf{x}\right)\propto\mathcal{CN}\left(\mathbf{x};\buildrel\lower 3.0pt\hbox{\scriptscriptstyle\rightarrow}\over{\mathbf{x}}_{c},\eta^{-1}_{c}\mathbf{I}\right) with $\eta_{c}=\omega_{c}-\tau_{c}$ and \buildrel\lower 3.0pt\hbox{\scriptscriptstyle\rightarrow}\over{\mathbf{x}}_{c}=\left(\omega_{c}\hat{\mathbf{x}}_{c}-\tau_{c}\boldsymbol{\gamma}_{c}\right)/\eta_{c}, which yields (8a) and (8b).

Subsequently, we select the factor $f_{0}$. In accordance with Step (2), we set the message ${\mu_{\mathbf{x}\to{f_{0}}}}\left(\mathbf{x}\right)$ to be the product of messages passing from factors $\{f_{1},f_{2},\ldots,f_{C}\}$ to $\mathbf{x}$, then

[TABLE]

Then, we have ${\mu_{\mathbf{x}\to{f_{0}}}}\left(\mathbf{x}\right)\propto\mathcal{CN}\left(\mathbf{x};\boldsymbol{\gamma}_{0},\tau^{-1}_{0}\mathbf{I}\right)$ with $\tau_{0}=\sum_{c\in\mathcal{C}}\eta_{c}$ and \boldsymbol{\gamma}_{0}=\sum_{c\in\mathcal{C}}\eta_{c}\buildrel\lower 3.0pt\hbox{\scriptscriptstyle\rightarrow}\over{\mathbf{x}}_{c}/\tau_{0}, which yields (9a) and (9b). Following Step (3), we set the approximate belief on $f_{0}$ as $b_{f_{0}}(\mathbf{x})\propto\mathcal{CN}\left(\mathbf{x};\hat{\mathbf{x}}_{0},\omega^{-1}_{0}\mathbf{I}\right)$ with $\omega_{0}=\sum_{k\in\mathcal{K}}v_{0,k}/K$, where $\hat{x}_{0,k}$ and $v_{0,k}$ denote the MMSE estimate and its corresponding MSE of $x_{k}$ from its AWGN observation modeled by (11), whose explicit expression are given by (12). Similar to (26), we set the message ${\mu_{{f_{0}}\to\mathbf{x}}}\left(\mathbf{x}\right)$ as

[TABLE]

Finally, we denote $\eta_{0}=\omega_{0}-\tau_{0}$ and \buildrel\lower 3.0pt\hbox{\scriptscriptstyle\leftarrow}\over{\mathbf{x}}_{0}=\left(\omega_{0}\hat{\mathbf{x}}_{0}-\tau_{0}\boldsymbol{\gamma}_{0}\right)/\eta_{0}, and the current iteration is completed. Afterwards, the next iteration begins, and the factors $\{f_{1},\ldots,f_{C}\}$ are selected again. From Step (2) and similar to (27), we compute the message ${\mu_{\mathbf{x}\to{f_{c}}}}(\mathbf{x})$ by

[TABLE]

Then, we can simply the expression of $\tau_{c}$ as

[TABLE]

which is identical to (6a). In the similar manner, we can simplify $\boldsymbol{\gamma}_{c}$ as (6b). Repeating the above message update procedure yields Algorithm 1.

Appendix B Proof of Proposition 1

We start by defining the following auxiliary functions:

[TABLE]

Then, we can rewrite the evolution equations in (13a) and (13b) w.r.t. $\{\Phi_{c}\}_{c=0}^{C}$ and $\Psi$ as follows can be re written

[TABLE]

where (30b) follows from the manipulation of (13b) shown as

[TABLE]

We learn from the appendices in [24] that $\Phi_{c}$’s and $\Psi$ are monotonically decreasing for nonnegative $\nu$ and $\rho$ and that

[TABLE]

As $\Psi\left(\rho\right)$ is monotonically decreasing, we conclude that for $\rho>0$,

[TABLE]

and evidently, $\Psi\left(\rho\right)>0$. We also notice that

[TABLE]

recalling that $\lambda_{i,c}\geq 0$, which follows directly that $\Phi_{c}\left(\nu\right)\geq 0$ for $\nu\geq 0$ and $\forall c\in\mathcal{C}$.

On the basis of the above results, we prove the monotonicity of the sequences $\{\rho^{t}\}$ and $\{\nu^{t}_{c}\}_{c=1}^{C}$ by induction. The monotonicity of $\{\nu^{t}_{c}\}$ can be proved by induction. Initially, we have

[TABLE]

where the inequality follows from (31) and the non-negativity of $\Phi_{c}$’s. Then, we proceed by supposing that $\nu^{t}_{c}<\nu^{t-1}_{c}$. As $\{\Phi_{c}\}_{c=1}^{C}$ are monotonically decreasing, we have $\Phi_{c}\left(\nu^{t}_{c}\right)>\Phi_{c}\left(\nu^{t-1}_{c}\right)$, which follows directly that

[TABLE]

Then from the monotonicity of $\Psi$, we have $\Psi\left(\rho^{t}\right)>\Psi\left(\rho^{t-1}\right)$. Integrating the fact that $\Phi_{c}\left(\nu^{t}_{c}\right)>\Phi_{c}\left(\nu^{t-1}_{c}\right)$, we find that

[TABLE]

Hence, $\{\nu^{t}_{c}\}_{c=1}^{C}$ are monotonically decreasing sequences. In addition, the monotonicity of $\{\rho^{t}\}$ can directly follow from (32).

For their boundedness, the non-negativity of $\Phi_{c}$’s and $\Psi$ implies that the sequences $\{\rho^{t}\}$ and $\{\nu^{t}_{c}\}$’s have the lower bound [math]. Meanwhile, their monotonicity implies they have upper bound $\nu^{1}_{c}$ and $\rho^{0}$, respectively.

Appendix C Proof of Proposition 2

In this appendix, we prove Proposition 2 via optimality condition for the equality constrained problem [35, Sec. 3.1]. We define the Lagrangian function corresponding to the optimization problem in (18) as

[TABLE]

To show that $b_{0}(\mathbf{x}),b_{1}(\mathbf{x}),\ldots,b_{C}(\mathbf{x})$ and $q(\mathbf{x})$ in (17) are the stationary points of the problem (18), we show that there exists Lagrange multipliers ${\boldsymbol{\alpha}}_{0},{\boldsymbol{\alpha}}_{1},\ldots,{\boldsymbol{\alpha}}_{C},{\boldsymbol{\beta}}_{0},{\boldsymbol{\beta}}_{1},\ldots,{\boldsymbol{\beta}}_{C}$ and ${\boldsymbol{\mu}}=[\mu_{0},\mu_{1},\ldots,\mu_{C}]^{T}$ such that $b_{0}(\mathbf{x}),b_{1}(\mathbf{x}),\ldots,b_{C}(\mathbf{x})$ and $q(\mathbf{x})$ in (17) satisfy

[TABLE]

Let ${\boldsymbol{\alpha}}_{c}=-\Re(\tau_{c}\boldsymbol{\gamma}_{c})$, ${\boldsymbol{\beta}}_{c}=-\Im(\tau_{c}\boldsymbol{\gamma}_{c})$, and $\mu_{c}=K\tau_{c}$, we can write the Lagrangian function w.r.t. $b_{c}$ as follows

[TABLE]

where the constant terms are independent on $b_{c}$, and the last two equalities follow directly from the definition of KL divergence. Therefore, $b_{0}(\mathbf{x}),b_{1}(\mathbf{x}),\ldots,b_{C}(\mathbf{x})$ defined in (17) achieves the minimization in (33a).

Similarly, we can obtain the maximization in (33b)by rewriting the Lagrangian function w.r.t. $q$ as follows

[TABLE]

where $(*)$ follows from Lemma 2 that $\sum_{c=0}^{C}\tau_{c}=C\omega$ and $\sum_{c=0}^{C}\tau_{c}\boldsymbol{\gamma}_{c}=C\omega\hat{\mathbf{x}}$. Hence, the maximization in (33b) can be attained by $q(\mathbf{x})$ defined in (17).

In addition, the moment-matching constraints in (19) can be justified to satisfy by noting that $\hat{\mathbf{x}}$ and $K\omega^{-1}$ are just the mean and trace of covariance matrices of the densities $b_{0}(\mathbf{x}),b_{1}(\mathbf{x}),\ldots,b_{C}(\mathbf{x})$ from (7b), (7c), (10a), and (10b) and that $q(\mathbf{x})$ is Gaussian with mean $\mathbf{x}$ and covariance matrix $w^{-1}\mathbf{I}_{K}$.

Appendix D Proof of Proposition 3

From the definition of the Stieltjes transform [36, Sec. 2.2.1], when $K\to\infty$,

[TABLE]

Then, we have $\tau_{c}=-S^{-1}_{\sigma^{-2}\mathbf{H}_{c}^{H}\mathbf{H}_{c}}(\omega^{-1}_{c})$, where $S^{-1}_{\sigma^{-2}\mathbf{H}_{c}^{H}\mathbf{H}_{c}}(\cdot)$ denotes the inverse function of $S_{\sigma^{-2}\mathbf{H}_{c}^{H}\mathbf{H}_{c}}(\cdot)$. Hence, from (8a) and the definition of R-transform [36, Sec. 2.2.5], we have

[TABLE]

According to Lemma 2, for any fixed point, $R_{\sigma^{-2}\mathbf{H}_{c}^{H}\mathbf{H}_{c}}(-\omega^{-1}_{c})=R_{\sigma^{-2}\mathbf{H}_{c}^{H}\mathbf{H}_{c}}(-\omega^{-1})$ for $\forall c\in\mathcal{C}$. We use $\hat{\mathbf{h}}_{c,i}$ and $\hat{\mathbf{h}}_{i}$ to denote the $i$-th column of the matrix $\mathbf{H}^{H}_{c}$ and $\mathbf{H}^{H}$ respectively. Then, following the argument in [36, Sec. 2.4.2] and from (9a), we can prove that when $\mathbf{H}$ is a random matrix with i.i.d. elements,

[TABLE]

The above arguments are derived from the asymptotic freeness of the relevant matrices, which also requires $K\to\infty$. As such, we complete the proof of (20).