A Low-Complexity Equalizer for Massive MIMO Systems Based on Array Separability

—In this paper, we propose a low-complexity equalizer for multiple-input multiple-output systems with large receive antenna arrays. The computational complexity reduction is achieved by exploiting array separability on a geometric channel model. This property suggests a two-stage receive processing, consisting of (i) sub-array beamforming and (ii) low-dimension minimum mean square error (MMSE) equalization. Simulations indicate that the proposed method outperforms the classical MMSE ﬁlter in terms of complexity provided that the number of channel scatterers and the sub-arrays dimensions are not excessively large.


I. INTRODUCTION
Mobile communication has established itself as a ubiquitous technology in modern life. The number of mobile connected devices is expected to grow exponentially owing to the advent of new technologies such as wireless sensor networks and to the popularity of services such as instant messaging and audio/video streaming. The quality of service in this demanding scenario relies on higher data-rates, broader signal coverage, and increased reliability from the communications systems. Advances in the physical layer are thus necessary in order to satisfy these continuously increasing system requirements. Among the most promising novel technologies for 5th generation of mobile communications systems are massive multiple-input multiple-output (MIMO) and mmWave technologies [1], [2]. The former enables wireless systems to form highly directive beams that can maximize energy transfer to a desired user, while minimizing leakage towards interfering users. The latter provides very large bandwidth and antenna miniaturization thanks to millimeter range propagation.
Linear filtering turns out to be a possible solution to the receive processing problem for massive MIMO systems due to their conceptual simplicity. For instance, the minimum mean square error (MMSE) filter aims at maximizing the signal portion while attenuating the interference and noise components. The major drawback of this filter is the complexity required for the data covariance matrix inversion, which may become prohibitive in the massive MIMO case. In view of this, strategies to reduce the computational complexity of linear filters are of interest.
Several complexity reduction strategies have been proposed as alternatives to the closed-form MMSE filter, including: Gauss-Seidel optimization [3], Neumann series approximation [4], and the symmetric successive overrelaxation method [5]. An approach that exploits the array manifold separability has been presented as a solution to reduce the computational complexity in large array beamforming tasks [6]. Therein, the Khatri-Rao product factorization of the array manifold matrix is exploited to simplify the beamforming filters calculation. However, [6] is limited to the single-input multiple-output case without multipath propagation.
Contribution: In this paper, we assume array separability at the receiver side to simplify the equalizer filter design in a massive MIMO system with two-dimensional arrays. As will be shown later, complexity reduction is achieved by decoupling sub-array beamforming and MMSE equalization into two stages. To this end, we introduce multilinear algebra notation in Section II and consider the Khatri-Rao factorization of a geometric channel model from the array separability assumption in Section III. The MIMO system model is recast using multilinear algebra to provide a spatially decoupled model. We then propose an equalization scheme based on multilinear algebra product properties in Section IV. Simulation results are shown and discussed in Section V, and the paper is concluded in Section VI.
Notation: We denote scalars by lower and uppercase italic letters, vectors by lowercase bold letters, matrices by uppercase bold letters, and tensors as calligraphic uppercase letters. The notation [·] (n) stands for the n-mode unfolding of the argument tensor. Matrix transpose and conjugate-transpose are respectively denoted by (·) T and (·) H . Frobenius and 2 norms are written as · F and · 2 , respectively. Symbols × n , ⊗ and stand for the n-mode, Kronecker and Khatri-Rao products, respectively. The statistical expectation operator and the pseudo-inverse are represented by E[·] and (·) † , respectively. Matrix I N denotes the (N × N ) identity matrix, and I N the (N × N × N ) identity tensor, whose values assume 1 along its main diagonal and 0 elsewhere.

II. MULTILINEAR ALGEBRA
Multilinear algebra is essentially the algebraic framework for multilinear operators, also known as tensors. In this work, the term "tensor" refers to a multidimensional array structure. The reader is referred to [7] for more formal definitions. In the following, we present definitions and properties for 3dimensional tensors since the treated structures in this work are, at most, 3-dimensional.
Consider a third-order tensor T ∈ C N1×N2×N3 . We define respectively the elements of the 1-, 2-, and 3-mode unfolding matrices of T : as the "matricization" of T along its three modes. The 1-, 2-, and 3-mode products between tensor T and matrices U 1 ∈ C M1×N1 , U 2 ∈ C M2×N2 and U 3 ∈ C M3×N3 , respectively, are defined as One can simultaneously apply n-mode products along different modes. For example, consider Equation (1) is known as the Tucker model, and, in this case, T is referred to as the core tensor and U 1 , U 2 , U 3 as its factor matrices.
A special case of the Tucker model is the Parallel Factors (PARAFAC) model, which is defined as Here, the core tensor is diagonal and the factor matrices A ∈ C M1×R , B ∈ C M2×R , C ∈ C M3×R have the same number R of columns, which is also called the tensor rank. The unfolding matrices of (2) can shown to be written as [8]: III. SYSTEM MODEL Consider a MIMO system with a transmitter equipped with a uniform planar array (UPA) of size N h t elements in the horizontal plane times N v t elements in the vertical plane The receiver employs a UPA of size N r = N h r N v r , as depicted in Figure 1. The horizontal and vertical element spacing for the transmit and receive arrays are δ h t , δ v t , δ h r , and δ v r , respectively. The transmitter applies a filter F ∈ C Nt×B to transmit the data symbol vector s(n) The receiver observes the following signal vector: where H ∈ C Nr×Nt is the narrow-band block fading channel matrix between the transmitter and receiver, and b(n) ∈ C Nr×1 the zero mean additive noise component with We adopt a geometric channel model with L scatterers, each one contributing a single propagation path between the transmitter and receiver. These assumptions are valid for transmissions in the mmWave band [9], [10]. The channel matrix can then be written as where λ is the complex gain of the th path, with E[|λ 2 |] = 1. The departure and arrival azimuth angles of the th path are denoted by φ t and φ r , respectively, whereas its departure and arrival elevation angles are denoted by θ t and θ r . The factors g t (·) and g r (·) refer to the antenna gains of the transmitting and receiving arrays, respectively. The response of the antenna element in the Cartesian position (n h t , n v t ) of the transmit array with respect to a plane wave from direction (θ, φ) can be written as [11] Likewise, the element response in position (n h r , n v r ) of the receive array is with n h r ∈ {1, . . . , N h r } and n v r ∈ {1, . . . , N v r }. Using matrix notation, Equation (5) can be rewritten as We also assume N r ≥ N t and B ≤ min(N r , N t ). The channel matrix H is assumed to be full column-rank, implying L ≥ N t .
It is well-known that the antenna array response can be separated into its spatial domain contributions [11]. Assuming separability on the receive array, it follows that ) and a v r (θ, φ) denote the steering vectors corresponding to the linear sub-arrays along the y and z-axes, respectively. In view of this, the receive array manifold matrix can be rewritten as Let H e HF ∈ C Nr×B be the effective channel matrix. It can be expressed as We identify (6) as the transposed 3-mode unfolding (3) of the third-order PARAFAC tensor Equation (7) represents a trilinear channel model that decouples the receive horizontal and vertical spatial signatures from the channel fading and transmitter beamforming. Note that (7) could be easily extended to support higher-order separable arrays, as shown in [6]. Using tensor notation, the received signal (4) is rewritten as where X(n) ∈ C N h r ×N v r and B(n) ∈ C N h r ×N v r denote reshaping matrices of x(n) and b(n), respectively. In the next section, we propose a novel equalization strategy that exploits the multilinearity of the data model (8) to reduce the computational cost of the receive signal processing.

IV. PROPOSED EQUALIZATION METHOD
The MMSE equalizer is one of the most popular MIMO receiver filters since it attempts to jointly minimize interference and noise contamination [11]. Its coefficients are obtained by computing As the number of receiving antennas grows, the calculation of G o becomes more convoluted due to the large computational costs involved in the inversion of R xx .
In order to reduce the equalizer's implementation costs, we propose to exploit the multilinear channel structure (7) in the equalizer optimization. It is assumed that the receiver knows the sub-arrays manifold matrices A h r and A v r and. In practice, these parameters can be estimated using array processing techniques such as MUSIC and ESPRIT [12]. Also, both the proposed and MMSE equalizers assume perfect knowledge of the channel matrix H. The proposed filtering method is divided into two stages: (i) sub-array beamforming, and (ii) lowdimension MMSE equalization. It will be hereafter referred to as Two-Step Sequential (2-Seq) filter.
More specifically, the sub-array beamforming task consists of performing spatial zero-forcing filtering along the first and second signal dimensions. The beamformer output is given by [8] The subsequent step consists of equalizing the fading matrix U in (10). To this end, we consider the vectorization of Y(n): the colored noise component. It is worth noting that in order to form Θ, one needs to first compute U = (A † r H e ) T . The MMSE equalizer for y(n) is then given bỹ Qo argmin From Equations (8), (10), (11), and (12), it follows that 2-Seq filter can be concisely written as In scenarios with poor scattering propagation (e.g. mmWave channels), computational complexity can be reduced further by capitalizing on the sparsity of Θ.
The number L of scatterers has a major role on the filter's numerical stability and computational complexity. The subarray manifold matrices A h r and A v r need to be full columnrank in order to their pseudo-inverse exist in (10). Since the probability of these matrices becoming singular increases as L grows, scenarios with only a few scatterers are more attractive for the use of the 2-Seq filter. Besides, the colored noise covariance matrix Rbb may become singular as the number L of scatterers grow, making R yy a near-singular matrix.
The proposed two-step filtering strategy can still be applied to the Rayleigh channels even though they do not lead to the multilinear representation in (8). In this case, the beamforming output signal is given by

The low-dimensional MMSE equalizer is then as followŝ
Qo argmin Here the sub-array manifold matrices are formed according to some strategy. We opted for pointing the array beams at L linear spaced angles in order to reduce the method's computational complexity. Nevertheless, any classical beamforming strategy for Rayleigh channels could still be employed, resulting in an increased complexity. It is worth noting that Φ does not exhibit a sparsity pattern as Θ because the Rayleigh channel model does not admit the multilinear representation (7). In this scenario, the 2-Seq filter is expected to be less accurate than the classical MMSE filter as a result of the model mismatch. However, the computational complexity reduction can be leveraged in practical scenarios.
Computational complexity: The first filtering stage corresponds to O(N h 3 r + N v 3 r ) operations due to the pseudoinverses. The computation of (12) or (13) amounts to O(L 6 ) operations since the beamformer output signal lies in a L 2dimensional space. Therefore, the general computational complexity of the 2-Seq filter is The calculation of the 2-Seq and MMSE filters amounts to O(3λ 3 ) and O(λ 6 ) operations, respectively. In the following section, results of numerical simulations conducted to evaluate the computational performance of the proposed method will be presented and discussed.

V. SIMULATION RESULTS
In this section, we investigate the computational complexity and the symbol recovery performance of the 2-Seq receiver. We also assess the influence of the channel model on its performance. The classical MMSE filter (9) is used as a benchmark.
We consider two figures of merit in our simulations: the number of FLOPS necessary to compute the equalizer and the MSE = 1 R R r=1 1 K S r − GX 2 F , where R denotes the number of independent Monte Carlo runs, K the number of transmitted QPSK symbols, S r ∈ C B×K the symbol matrix in the rth experiment, X ∈ C Nr×K the corresponding received signals matrix, and G the considered equalizer. The number of FLOPS were counted using the Lightspeed toolbox [13]. In each experiment, the channel gains λ were generated from a zero mean and unit variance Gaussian random variable, the direction of departure and arrival azimuth and elevation angles were uniformly selected from [−80 • , 80 • ] and [10 • , 170 • ], respectively. The horizontal and vertical array inter-element spacing is half-wavelength. Truncated discrete Fourier transform matrices were used for precoding.
In Figure 2, we provide the MSE performance and the mean processing time for the MMSE and 2-Seq filters as functions of the number N r of antenna elements in the receive array for R = 1000 experiments and L ∈ {6, 9, 12} scatterers. The receive array consists of a UPA with dimensions N r = N × N, N ∈ {15, 20, 25, 30}. Since the N > L in all scenarios, the sub-array manifold matrices are tall and possibly full column-rank depending on the multipath arrival angles. The very few experiments where the 2-Seq filter failed due to numerical instability were discarded according to a MSE threshold. Figure 2 indicates that the symbol recovery quality of both filters improves as L and N r increase since there are more "copies" of the transmitted data. By contrast, 2-Seq's computational complexity varies little with N r , as expected from the discussion in Section IV. This result suggests that the proposed method becomes more efficient than the benchmark one as L decreases, which may be the case in scenarios with poor scattering propagation.
The influence of the noise power on the filtering performance is investigated in Figure 3 for the geometrical (5) and Rayleigh channel models, where the results were obtained from R = 1000 experiments. The inputs of the Rayleigh channel matrix were extracted from i.i.d. zero mean and unit variance Gaussian random variables for each experiment. We observe in the geometrical model scenario that there is no MSE degradation for the 2-Seq filter due to the system separability.
For the Rayleigh channel model, however, there is a MSE performance gap. In this case, H does not admit the separable representation (7), and the signal model (8) becomes inaccurate. Nevertheless, the proposed method is still employed in this experiment by forming the beamformers A h r and A v r pointed at L = 6 linearly spaced directions and afterwards performing the low-dimension MMSE equalization (13). The performance gap is therefore due to the low-dimensionality of the MMSE filter compared to the Rayleigh channel rank. The number of FLOPS is practically the same for both channel scenarios, since the same value for N r and L were used. The small difference comes from the computation of Θ and Φ.
The absolute value of the colored noise covariance matrix inputs are depicted in Figure 4 for L = 6 and L = 12. We observe that as L grows, Rbb becomes sparser and rank deficient and, consequently, Θ approaches singularity. Besides the low-complexity, numerical stability is therefore another reason to employ the proposed method in communication channels with only a few scatterers.

VI. CONCLUSION
We presented a low-complexity equalization scheme for MIMO systems with massive receive antenna array. The proposed filter is based on the assumed receive array separability reflecting a multilinear structure on the geometrical channel and data models. Such algebraic property is exploited to decouple spatial filtering from MMSE equalization, which enables computational cost reduction. Simulation results confirm the reduced complexity of the proposed method provided the number of channel scatterers is small. Furthermore, results indicate that the proposed filter can still be employed to nonseparable channels resulting in computational cost reduction with symbol recovery degradation. Full transceiver optimization considering separability at both transmitter and receiver ends with interfering users, and channel estimation methods based on the presented model are outlined as future work.