{'fraction_non_alphanumeric': 0.03618280030866844, 'fraction_numerical': 0.009002829460687645, 'mean_word_length': 4.452664796633941, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Enforcing Statistical Orthogonality in Massive MIMO Systems via Covariance Shaping

This paper tackles the problem of downlink data transmission in massive multiple-input multiple-output (MIMO) systems where user equipments (UEs) exhibit high spatial correlation and channel estimation is limited by strong pilot contamination. Signal subspace separation among UEs is, in fact, rarely realized in practice and is generally beyond the control of the network designer (as it is dictated by the physical scattering environment). In this context, we propose a novel statistical beamforming technique, referred to as MIMO covariance shaping, that exploits multiple antennas at the UEs and leverages the realistic non-Kronecker structure of massive MIMO channels to target a suitable shaping of the channel statistics performed at the UE-side. To optimize the covariance shaping strategies, we propose a low-complexity block coordinate descent algorithm that is proved to converge to a limit point of the original nonconvex problem. For the two-UE case, this is shown to converge to a stationary point of the original problem. Numerical results illustrate the sum-rate performance gains of the proposed method with respect to spatial multiplexing in scenarios where the spatial selectivity of the base station is not sufficient to separate closely spaced UEs.

Abstract-This paper tackles the problem of downlink data transmission in massive multiple-input multiple-output (MIMO) systems where user equipments (UEs) exhibit high spatial correlation and channel estimation is limited by strong pilot contamination. Signal subspace separation among UEs is, in fact, rarely realized in practice and is generally beyond the control of the network designer (as it is dictated by the physical scattering environment). In this context, we propose a novel statistical beamforming technique, referred to as MIMO covariance shaping, that exploits multiple antennas at the UEs and leverages the realistic non-Kronecker structure of massive MIMO channels to target a suitable shaping of the channel statistics performed at the UE-side. To optimize the covariance shaping strategies, we propose a low-complexity block coordinate descent algorithm that is proved to converge to a limit point of the original nonconvex problem. For the two-UE case, this is shown to converge to a stationary point of the original problem. Numerical results illustrate the sum-rate performance gains of the proposed method with respect to spatial multiplexing in scenarios where the spatial selectivity of the base station is not sufficient to separate closely spaced UEs.

I. INTRODUCTION
Massive multiple-input multiple-output (MIMO) is a multiantenna technology that has great potential to boost the spectral efficiency (SE) of cellular networks by means of highly directional beamforming and spatial multiplexing of many user equipments (UEs) in the same time-frequency resources. It thus plays a pivotal role in current 5G New Radio (NR) implementations [3]- [5] and is expected to maintain this prominence in future wireless generations [6]- [10]. The benefits of massive MIMO can be ascribed to the large number of antennas available at the base station (BS), which we denote by M . In this context, it is shown in [11] that the achievable SE of downlink/uplink massive MIMO systems is unbounded as M grows large and when minimum mean squared error (MMSE) precoding/combining is adopted at the BS, which can asymptotically remove any interference. However, in presence of a large number of UEs and finite BS antennas, the aforementioned approach might still be limited by interference and must rely on accurate instantaneous channel state information (CSI). In crowded scenarios, such as outdoor events, transport hubs, and stadiums, the channels of closely spaced UEs exhibit high spatial correlation, which hinders the capability of the BS to separate such UEs during both the channel estimation phase and the data transmission phase. Moreover, when the channel coherence time is limited, non-orthogonal pilots might be used across the UEs during the channel estimation phase, which results in strong pilot contamination [12], [13].
To overcome these issues and facilitate the operations in the massive MIMO regime, several works have proposed to leverage statistical CSI [14]- [25]. This consists mainly in the channel covariance matrix of each UE, which is essentially dictated by the angle spread spanned by the multipath propagation of the signals impinging on the antenna array. This angle spread is often bounded due to the high spatial resolution of the massive array compared with the limited scattering. As a result, the channel covariance matrices in massive MIMO tend to be low-rank and dominated by few major propagation directions [11], [17], [18]. This particular property can be exploited for several applications such as reducing the feedback overhead in the channel estimation phase [14]- [16] and mitigating the interference in the downlink data transmission phase when the UEs exhibit non-overlapping or orthogonal signal subspaces [17], [18]. For example, statistical CSI can be used to precode signals such that their chosen propagation paths do not interfere in average [19], [20]. In this regard, a robust precoding/decoding design based on the average mean squared error (MSE) matrix is derived in [19] under different CSI conditions, whereas [20] provides a lower bound on the ergodic sum rate for the two-UE setting when the UEs are equipped with a single antenna.
There is also a large body of literature focusing on the hybrid precoding problem, where the beamforming applied at the BS is factorized into an inner and an outer precoding matrix [14]- [16], [21], [22]. Here, the former is based on instantaneous CSI while the latter depends only on the secondorder channel statistics. This results in reduced pilot length required for the instantaneous effective channel estimation and lower computational complexity associated with the inner precoding design. Moreover, in the context of frequency-division duplex (FDD) systems, such a hybrid precoding method can be used to reduce the size of the resulting effective channels arXiv:2106.07952v2 [cs.IT] 30 Mar 2022 by exploiting the near-orthogonality of the angle spreads of different UE groups [14]. In [15], this approach is shown to incur a vanishing loss compared with the case of full CSI when the number of BS antennas grows large, while the same result is obtained in [16] for the case of finite number of BS antennas and suitable UE scheduling. In time-division duplex (TDD) systems, existing works have pointed out how the exploitation of statistical CSI enables efficient pilot reuse across the UEs [23]- [25]. In [23], pilot reuse for massive MIMO transmissions over spatially correlated Rayleigh fading channels is proposed when the angle spreads of the UEs with the same pilots are non-overlapping. In [24], the second-order channel statistics are used to precode the pilots, thus reducing the variance of the channel estimation error by a factor that is proportional to the number of UE antennas. A robust channel estimation method dealing with scenarios in which the angle spreads of the desired and interfering channels are overlapping is proposed in [25] and exploits channel separability in the power domain [26].
Most of the above works rely on specific structures of the channel covariance matrices of the UEs in terms of rank and degree of separation among signal subspaces. Such properties are determined by the locations of the UEs combined with the surrounding scattering environments, where both these factors are generally beyond the control of the network designer. Indeed, in many practical scenarios, the UEs exhibit high spatial correlation, e.g., when they are not sufficiently far apart. Hence, the condition of non-overlapping or orthogonal signal subspaces is rarely satisfied in practice [11]. Note that, in this setting, the number of BS antennas needed to effectively separate interfering UEs with MMSE precoding/combining increases with their spatial correlation, i.e., as the UEs get closer to each other, and may become impractical in crowded scenarios. Moreover, existing works are often based on the assumption of Kronecker channel model owing to its simple analytical formulation. However, such a model has been shown to be an over-simplification of the true nature of MIMO channels [27]. 1

A. Contributions
Building on the fact that most future UEs are envisioned to be equipped with a small-to-moderate number of antennas, this paper proposes a novel statistical beamforming technique at the UEs-side, referred to as MIMO covariance shaping. 2 While existing methods assume statistical orthogonality among the UEs as a property given by the physical scattering environment and the Kronecker channel model due to its analytical tractability, MIMO covariance shaping aims at modifying the channel statistics of the UEs with highly overlapping channel covariance matrices in order to enforce signal subspace separation in scenarios where the spatial selectivity of the BS is not sufficient to separate such UEs. The proposed approach consists in preemptively applying a statistical beamforming at the UE-side during both the uplink pilot-aided channel estimation phase and the downlink data transmission phase, and relies uniquely on statistical CSI. In this context, each UE employs its antennas to excite only a subset of all the possible propagation directions towards the BS such that the spatial correlation among the interfering UEs is minimized while preserving enough useful power for effective data transmission. Hence, MIMO covariance shaping is suitable for both pilot decontamination in TDD systems and statistical precoding/combining. Remarkably, the proposed method exploits the realistic non-Kronecker channel structure, which allows to suitably alter the channel statistics perceived at the BS by designing the transceiver at the UE-side. Therefore, it has the unique advantage of turning the generally inconvenient non-Kronecker nature of massive MIMO channels into a benefit. Numerical results show the sum-rate performance gains with respect to a reference scheme employing multiple antennas at the UE for spatial multiplexing in scenarios where the spatial selectivity of the BS is not sufficient to separate UEs with highly overlapping channel covariance matrices and the channel estimation is limited by strong pilot contamination.
The MIMO covariance shaping framework was initially proposed in our prior works [1], [2] and recently used for the minimization of the outage probability in [30]. This work extends the metric proposed in [1] to all pairs of interfering UEs to design a statistical receive beamforming at the UE along with a statistical precoding at the BS. Differently from [30], we argue that MIMO covariance shaping must be adopted for both pilot decontamination during the uplink pilot-aided channel estimation phase and combining during the downlink data transmission phase. By doing so, the BS is able to acquire accurate instantaneous estimates of the effective channels, which are then used for efficient downlink precoding and result in a great enhancement of the network performance.
The contributions of this paper are summarized as follows.
• We present the novel concept of MIMO covariance shaping, which aims at designing a suitable shaping of the channel covariance matrices of the UEs to enforce a full or partial separation of their signal subspaces that would be otherwise highly overlapping. • We point out that the exploitation of the non-Kronecker nature of massive MIMO channels is crucial to suitably alter the channel statistics perceived at the BS by acting at the UE-side. • We derive a tractable expression of the ergodic achievable sum rate under MIMO covariance shaping with the objective of characterizing the impact of the proposed framework on the system performance. • We optimize the covariance shaping strategies by minimizing the variance of the inter-UE interference (as a metric to measure the spatial correlation) among the interfering UEs. To this end, we propose a low-complexity block coordinate descent algorithm that is proved to converge to a limit point of the original nonconvex problem. For the two-UE case, this is shown to converge to a stationary point of the original problem. • We provide numerical results characterizing several scenarios where MIMO covariance shaping outperforms spatial multiplexing. In particular, a superior sum-rate performance is achieved when the spatial selectivity of the BS is not sufficient to separate UEs exhibiting high spatial correlation and the channel estimation is limited by strong pilot contamination. Outline. The rest of the paper is structured as follows. Section II introduces the system model for a general multi-UE MIMO downlink system. Section III presents the concept of MIMO covariance shaping. Section IV proposes efficient methods to optimize the covariance shaping vectors. Section V presents numerical results to evaluate the performance of the proposed framework. Finally, Section VI concludes the paper.
Notation. Lowercase and uppercase boldface letters denote vectors and matrices, respectively, whereas (·) T , (·) H , and (·) * are the transpose, Hermitian transpose, and conjugate operators, respectively. · and · F represent the Euclidean norm for vectors and the Frobenius norm for matrices, respectively, whereas E[·] and V[·] are the expectation and variance operators, respectively. I A denotes the A-dimensional identity matrix and 0 represents the zero vector or matrix with proper dimension. tr(·) and vec(·) are the trace and vectorization operators, respectively. (a mn ) M m,n=1 denotes the matrix of size M × M whose (m, n)th element is given by a mn , whereas [a 1 , . . . , a A ] represents horizontal concatenation. We use ⊗ to denote the Kronecker product and u min [A] to denote the eigenvector corresponding to the minimum eigenvalue of matrix A, with u min (A) = 1. 1 A denotes the indicator function, which is equal to 1 if condition A is satisfied and to 0 otherwise, and · is the ceiling operator. Lastly, CN (0, A) is the circularly symmetric complex Gaussian distribution with zero mean and covariance matrix A.

II. SYSTEM MODEL
In this section, we introduce the channel model for a general multi-UE MIMO downlink system. Then, we describe the uplink pilot-aided channel estimation assuming a TDD setting and channel reciprocity between the uplink and the downlink. Finally, we discuss the system model for the downlink data transmission.

A. Channel Model
Consider a multi-UE MIMO system where a BS equipped with M antennas serves K UEs with N antennas each in the where h k,m ∈ C N ×1 and g k,n ∈ C 1×M are the channel vectors between the mth BS antenna and UE k and between the BS and the nth antenna of UE k, respectively. We assume a correlated Rayleigh fading channel model where the entries of Here, the channel covariance matrix Σ k ∈ C N M ×N M has the following general structure: 22 . . .
where each block Σ k,mn E[h k,m h H k,n ] ∈ C N ×N represents the cross-covariance matrix between the mth and nth columns of H k . Lastly, we define the covariance matrix seen at UE k as R k E[H k H H k ] ∈ C N ×N and the covariance matrix relative to UE k seen at the BS as T k E[H H k H k ] ∈ C M ×M , respectively. Observe that, in the case of downlink data transmission, R k and T k represent the receive and transmit covariance matrices, respectively.

B. Uplink Pilot-Aided Channel Estimation
Assuming a TDD setting and channel reciprocity between the uplink and the downlink (see, e.g., [32, Ch. 1.3.5]), the channel matrices {H k } K k=1 are estimated at the BS using antenna-specific uplink pilots, such that N pilot vectors per UE are required. Let S p {k : UE k has pilot P p } be the set of UEs sharing the same pilot matrix P p ∈ C N ×τ , with p = 1, . . . , P . Here, P < K denotes the number of orthogonal pilot matrices and τ represents the pilot length. The orthogonal pilot matrices satisfy {P p P H p = τ N I N } P p=1 , i.e., N orthogonal pilot vectors are assigned to each UE, and {P p P H q = 0} p =q . Note that these conditions imply τ ≥ P N . Furthermore, Y ∈ C M ×τ denotes the receive signal at the BS during the uplink pilot-aided channel estimation phase, which is given by where ρ UE is the transmit power at the UEs and Z ∈ C M ×τ is the noise term at the BS with elements independently distributed as CN (0, σ 2 BS ). Let Φ k,nn E[g T k,n g * k,n ] ∈ C M ×M denote the covariance matrix of g k,n . Then, the MMSE estimate of g k,n , with k ∈ S p , is given by (see, e.g., [32, where e n ∈ R N ×1 is the nth column of I N , Q k,nn Φ k,nn + j∈Sp\{k} Φ j,nn + N τ UE I M ∈ C M ×M is the normalized covariance matrix of the receive signal after correlation with P H p e n , and UE ρUE σ 2

BS
. Note that (3) can be seen as a superposition of channels estimated using the same pilot P H p e n , which cannot be distinguished by the BS: this phenomenon is referred to as pilot contamination [12].

C. Downlink Data Transmission
In the following, we focus on the downlink data transmission and assume that the BS transmits L k ≤ N independent symbols to UE k. We denote the transmit symbol vector for UE k by s k ∈ C L k ×1 , with E[s k s H k ] = I L k , and introduce the multi-UE transmit symbol vector s Before the transmission, the BS precodes s using the multi-UE precoding matrix is the precoding matrix corresponding to s k and W satisfies the power constrain W 2 F = 1. The receive signal at UE k is thus given by where ρ BS is the transmit power at the BS and z k ∼ Assuming perfect CSI at the UEs, the resulting sum rate is given by with BS ρBS σ 2

III. COVARIANCE SHAPING AT THE UE-SIDE
Although MMSE precoding/combining can asymptotically remove any interference as M grows large [11], this no longer holds in the presence of a large number of UEs and finite BS antennas. In particular, the number of BS antennas needed to effectively separate interfering UEs with MMSE precoding/combining increases with their spatial correlation, i.e., as the UEs get closer to each other, and may become impractical in crowded scenarios. In this setting, the BS can spatially separate signals corresponding to different UEs and mitigate or eliminate pilot contamination if their channel covariance matrices lie on orthogonal supports, i.e., if Σ k Σ j = 0 for a given pair of UEs k and j (see, e.g., [11]). However, the degree of statistical orthogonality among the UEs is determined by their locations combined with the surrounding scattering environments and both these factors are generally beyond the control of network designers. Hence, signal subspace separation among the UEs rarely occurs in practice.
In this context, we propose a novel method relying uniquely on statistical CSI and referred to as MIMO covariance shaping (or simply covariance shaping), which is applied at the UEside to enforce the aforementioned signal subspace separation in scenarios where the spatial selectivity of the BS is not sufficient to separate the UEs. While covariance shaping is conceived in such a way that it is agnostic to the operating frequency, we envision its application especially in sub-6 GHz frequency bands, where the availability of BS antennas may be limited and further exacerbates the lack of sufficient spatial selectivity. According to covariance shaping, the UEs preemptively apply a statistical beamforming vector, different for each UE, that aims at spatially separating their transmissions, thus drastically reducing both pilot contamination and interference. Here, the original MIMO channel of each UE is transformed into an effective multiple-input single-output (MISO) channel by combining the transmit/receive signal with the corresponding statistical beamforming vector. To this end, the knowledge of the channel covariance matrices of the UEs is required. Although the time scale at which such information must be acquired and used to compute the covariance shaping vectors increases with the mobility of the UEs (or, more generally, as the channel coherence time reduces), this still happens much less frequently than the instantaneous uplink pilot-aided channel estimation. Remarkably, the proposed method exploits the realistic non-Kronecker structure of massive MIMO channels that allows to suitably alter the channel statistics perceived at the BS by acting at the UE-side, thus turning a generally inconvenient model into a benefit.
Let v k ∈ C N ×1 denote the statistical beamforming vector preemptively applied at UE k: in the rest of the paper, we refer to v k as covariance shaping vector. 3 Hence, the original MIMO channel H k between the BS and each UE k is with Σ k introduced in (1) and where E[ ḡ k 2 ] = tr(Φ k ). In the rest of this section, we describe how the two phases of uplink pilot-aided channel estimation and downlink data transmission are modified under covariance shaping and we provide a tractable expression of the resulting ergodic achievable sum rate. The optimization of the covariance shaping vectors is discussed in Section IV. For ease of exposition, we focus on the case where covariance shaping is applied to all the UEs in the systems. Nonetheless, the proposed method may be applied to one or more subsets of closely spaced UEs, e.g., those who cannot be effectively separated with MMSE precoding/combining. In this context, UEs adopting covariance shaping and others employing the multiple antennas for spatial multiplexing can be served simultaneously by the BS in a transparent manner.

A. Uplink Pilot-Aided Channel Estimation
To estimate the effective channels {ḡ k } K k=1 resulting from covariance shaping, the BS assigns the same pilot vector p p ∈ C 1×τ to all UEs k ∈ S p , with { p p 2 = τ } P p=1 and {p p p H q = 0} p =q . The receive signal at the BS during the uplink pilotaided channel estimation phase, which we denote byȲ ∈ C M ×τ , is given by (cf. (2)) Then, the MMSE estimate ofḡ k , with k ∈ S p , is given by where Q k Φ k + j∈Sp\{k}Φ j + 1 τ UE I M ∈ C M ×M is the normalized covariance matrix of the receive signal after correlation with p H p . Note that the estimation of the effective channels only requires one pilot vector per UE. Hence, the application of covariance shaping can potentially reduce the pilot length with respect to the estimation of the channel matrices {H k } K k=1 described in Section II-B, which requires N pilot vectors per UE.

B. Downlink Data Transmission
Focusing on the downlink data transmission, the BS now transmits only one symbol s k ∈ C to each UE k, i.e., {L k = 1} K k=1 . Hence, we have s = [s 1 , . . . , s K ] T ∈ C K×1 and the multi-UE precoding matrix becomes W = [w 1 , . . . , w K ] ∈ C M ×K , where w k ∈ C M ×1 is the precoding vector corresponding to s k . The receive signal at UE k is thus given by (cf. (4)) and, assuming perfect CSI at the UEs, the resulting sum rate is given by (cf. (5))

C. Ergodic Achievable Sum Rate
We now analyze the sum rate as a function of the covariance shaping vectors with the objective of characterizing the impact of the proposed framework on the system performance. In particular, we derive a tractable expression of the ergodic achievable sum rate by assuming that each UE k does not know the effective scalar channelḡ k w k = v H k H k w k instantaneously but only its expected value. Building on [32,Thm. 4.6], we can express the receive signal in (10) as where the effective noise term accounts for the lack of instantaneous CSI in addition to the interference and the noise at UE k. A lower bound on the signal-to-interferenceplus-noise ratio (SINR) achievable by UE k can be obtained by treating z k as Gaussian noise, which yields the effective SINR The expression in (13) can be further simplified by considering the case where the BS adopts maximum ratio transmission (MRT) precoding. In this setting, the multi-UE precoding matrix can be written as whereH [ḡ T 1 , . . . ,ḡ T K ] T ∈ C K×M andĤ [ĝ T 1 , . . . ,ĝ T K ] T ∈ C K×M contain the effective channels and their MMSE estimates defined in (9), respectively. Hence, when MRT precoding is adopted at the BS, the effective SINR of UE k in (13), with k ∈ S p , reduces to This can be rewritten in closed form as in (16), shown at the top of this page, by substituting the following terms: . We refer to Appendix A for the detailed derivations. Note that, in the case of perfect channel estimation, i.e., when UE → ∞ and all the UEs have orthogonal pilots, we have {Φ k Q −1 kΦ k =Φ k } K k=1 and 1 j∈Sp = 0 j =k . In this context, the effective SINR in (16) simplifies as Finally, Proposition 1 presents an ergodic achievable sum rate with MRT in the two cases of perfect and imperfect channel estimation. The proof follows from [32,Thm. 4.6] and is thus omitted.
Proposition 1. Assume that the BS adopts MRT precoding. Then, an ergodic achievable sum rate is given bȳ with γ k defined in (16) and in (17) for imperfect and perfect channel estimation, respectively. Furthermore, the following limit holds as M → ∞: withR defined in (11).

IV. COVARIANCE SHAPING OPTIMIZATION
In this section, we address the sum rate maximization through a proper design of the covariance shaping vectors at the UEs. To this end, we consider the variance of the inter-UE interference as a metric to measure the spatial correlation or, in other words, the degree of statistical orthogonality, between two interfering UEs. In the following, we first consider the simple case of K = 2 and then extend the resulting analysis to the general case of K ≥ 2.
Let us define the inter-UE interference between UEs k and j after applying covariance shaping as where the notationḡ k (v k ) makes explicit the dependence of the effective channelḡ k on the corresponding covariance shaping vector v k . The effective channels of UEs k and j yield asymptotic favorable propagation if they satisfy [ For a practical number of BS antennas M , a meaningful performance metric is the variance of Ω(v k , v j ), expressed as where the notationΦ k (v k ) makes explicit the dependence of the effective channel covariance matrixΦ k on the corresponding covariance shaping vector v k . Note that the terms in the numerator and denominator of (23) can be written as where we recall that Σ k,mn represents the cross-covariance matrix between the mth and nth columns of the original channel H k , i.e., before applying covariance shaping, as defined in (1). Observe that δ(v k , v j ) = 0 impliesΦ kΦj = 0, i.e.,Φ k andΦ j lie on orthogonal supports [11]. Indeed, considering the eigenvalue decomposition of the effective channel covariance matrices, which may be written as {Φ i = U i Λ i U H i } i∈{k,j} , the condition of statistical orthogonality is satisfied if U k = U j and tr(Λ k Λ j ) = 0, which implies the rank deficiency of bothΦ k andΦ j . Clearly, in the general case, this imposes M 2 conditions whereas only 2N variables can be adjusted: this means that the resulting system of equations can be solved when N ≥ M 2 2 , which is generally not verified in practice since M KN in massive MIMO scenarios. Hence, since full signal subspace separation can hardly be achieved by simple transceiver design at the UE-side, it is of interest to minimize the spatial correlation between each pair of interfering UEs.

A. Two-UE Case
In the two-UE case, i.e., for K = 2, the covariance shaping vectors of UEs k and j are computed by solving the optimization problem Algorithm 1 (Covariance shaping: Alternating optimization algorithm) j , Σ k , Σ j , and the accuracy . Fix n = 1.  (29). (23). Although problem (P1) is not convex in either v k or v j , a suboptimal solution can be efficiently obtained via alternating optimization, as suggested in [1]. Let us introduce the definition The optimal covariance shaping vector of UE k for a given v j , denoted by v k , is obtained as Since (28)  η j,mn (v j )Σ k,mn (29) and the optimal covariance shaping vector of UE j for a given v k is obtained in a similar way. Hence, problem (P1) is solved by alternating the optimization between v k and v j until a predetermined convergence criterion is satisfied, e.g., until the difference between the values of the objective in consecutive iterations is sufficiently small. This scheme is formalized in Algorithm 1, whose convergence properties are characterized in Proposition 2.
Proposition 2. The alternating optimization algorithm described in Algorithm 1 converges to a stationary point of problem (P1).
Proof: Observe that the objective function in (28) does not depend on the scaling of v k and, if we relax the nonconvex constraint v k 2 = 1 as v k 2 ≤ 1, (29) remains a solution of (28). Then, according to [34, Corollary 2], every limit point of the sequence generated by the alternating optimization algorithm applied to problem (P1) with the relaxed constraints is a stationary point of the original problem.
Algorithm 1 can be implemented in a centralized manner at the BS and the covariance shaping vectors are fed back to the corresponding UEs. Alternatively, it can be implemented  in a distributed fashion at the UEs. In this case, each UE can compute its covariance shaping vector without any information exchange with the other UE provided that the channel statistics of the latter are known, the order of update is fixed, and the same initial points are used. 4 In fact, under these conditions, the entire alternating optimization procedure can be carried out locally and independently at each UE. Note that the same considerations hold for Algorithm 2 presented in the next section for the multi-UE case, where we also discuss the computational complexity. The result of the covariance shaping optimization heavily depends on the physical scattering environment. Consider the non-line-of-sight (NLoS) scenario in Fig. 1(a), where there is no line-of-sight (LoS) path between the BS and the UEs, the UEs are equipped with uniform linear arrays (ULAs), and the inter-UE distance is d = 4 m. In this case, as shown in Fig. 1(b), the covariance shaping vectors tend to focus their power along reflected paths that interfere with each other as little as possible while also carrying sufficient channel power, which results in a nearly interference-free 4 Note that v   transmission/reception. On the other hand, when a LoS path exists, as in the LoS scenario in Fig. 2(a), this generally carries more channel power than any other path. In this case, as shown in Fig. 2(b), the propagation directions selected by the covariance shaping vectors tend to partially capture the LoS path while also focusing some of their power along separated reflected paths in order to achieve some degree of statistical orthogonality between the UEs. Remarkably, in both Fig. 1(b) and Fig. 2(b), the result of Algorithm 1 is very close to the optimal solution of problem (P1) obtained with exhaustive search within a set of 10 3 uniformly distributed vectors over the N -dimensional unit sphere for each UE, where the former is characterized by negligible complexity with respect to the latter.
Kronecker Channel Model. Let us consider the particular case where each channel H k follows the Kronecker channel model [28]. In this setting, we have H k = R  CN (0, I N M ). Accordingly, the channel covariance matrix in (1) can be expressed as Σ k = T T k ⊗ R k with block elements given by Σ k,mn = T * k,mn R k , where T k,mn denotes the (m, n)th element of T k . Hence, from (23), we have It is straightforward to observe that, in this case, δ(v k , v j ) is independent of v k and v j . Hence, under the Kronecker channel model, it is not possible to alter the channel statistics perceived at one end of the communication link by designing the transceiver at the other end: as a consequence, no meaningful effective channel separation can be performed. This is in accordance with the properties of the Kronecker channel model, whereby the transmit and receive covariance matrices are independent. In this context, the signal subspace separation is exclusively determined by the physical scattering environment and can only be achieved when tr(T k T j ) = 0, which is rarely satisfied in practice [27], [28].

B. Multi-UE Case
In the general case, i.e., for K ≥ 2, the covariance shaping vectors for each UE are computed by solving the optimization problem (23). Although problem (P2) is not convex in any of the optimization variables {v k } K k=1 , a suboptimal solution can be efficiently obtained via block coordinate descent, which can be interpreted as an extension of the alternating optimization approach presented in the previous section for the two-UE case. The optimal covariance shaping vector of UE k for given {v j } j =k is obtained as with η j,mn (v j ) defined in (27). Like (28), (33) is in the form of generalized Rayleigh quotient and thus admits the solution and the optimal covariance shaping vectors of the other UEs j for given {v k } k =j are obtained in a similar way. Hence, problem (P2) is solved by optimizing the strategy of each UE given the strategies of the other UEs until a predetermined convergence criterion is satisfied. Furthermore, at each iteration i, the update v can be used to Algorithm 2 (Covariance shaping: Block coordinate descent algorithm) , α ∈ (0, 1], and the accuracy . (27), ∀j = k. S.2: Compute v k as in (34).
limit the variation of the covariance shaping vectors between consecutive iterations, where the step size α ∈ (0, 1] is chosen to strike the proper balance between convergence speed and accuracy (see, e.g., [10], [35]). The proposed scheme is formalized in Algorithm 2, whose convergence properties are characterized in Proposition 3.
Proposition 3. The block coordinate descent algorithm described in Algorithm 2 converges to a limit point of problem (P2).
Proof: At each iteration i of Algorithm 1, the covariance shaping vector of UE k results from solving (28), which admits the optimal solution in (29). Hence, the sequence δ(v j ) i converges to a finite non-negative value.
The computational complexity of Algorithm 2 is essentially dictated by step S.2, which corresponds to the computation of the optimal covariance shaping vector of UE k for given {v j } j =k . This is expressed in (34) and has a complexity O(N 2.37 ) due to the eigenvalue decomposition and inverse matrix operations [25]. Hence, the overall computational complexity of Algorithm 2 is O(IKN 2.37 ), where I is the number of iterations required for convergence. Note that both N and I are fairly modest since the UEs are equipped with a smallto-moderate number of antennas and the algorithm converges in very few iterations (i.e., less than 10 in the simulation scenarios considered in Section V). We can thus conclude that Algorithms 1 and 2 are remarkably efficient in terms of computation complexity.

V. NUMERICAL RESULTS AND DISCUSSION
In this section, we present numerical results to evaluate the gains achieved by the proposed covariance shaping framework.
To this end, we compare the following alternative transmission/reception schemes, where the first is based on covariance shaping and the second is considered as a reference. 5 (A) Covariance shaping. The multiple antennas at the UE are employed to implement covariance shaping and one data stream per UE is transmitted using MMSE precoding. Specifically, each UE k applies its covariance shaping vector v k , obtained with Algorithm 1 or Algorithm 2 (depending on the value of K), during both the uplink pilot-aided channel estimation phase and the downlink data transmission phase, as discussed in Section IV. The BS obtains the MMSE estimates of the resulting effective channels {ḡ k } K k=1 as in (9) based on the transmission of P < K orthogonal pilot vectors by the UEs. These estimates are used to compute the (M × K)-dimensional MMSE precoding matrix at the BS as in [11]. (B) Spatial multiplexing. The multiple antennas at the UE are employed for spatial multiplexing and N data streams per UE are transmitted using MMSE precoding. Specifically, the BS obtains the MMSE estimates of the channel matrices {H k } K k=1 as in (3) based on the transmission of P < K orthogonal pilot matrices by the UEs. These estimates are used to compute the (M × KN )-dimensional MMSE precoding matrix at the BS as in [11]. Some comments are in order.
• The reference scheme in (B) requires antenna-specific pilots for the estimation of the channel matrices {H k } K k=1 . In particular, N orthogonal pilot vectors are assigned to each UE to avoid any pilot contamination among the antennas of the same UE, which implies τ ≥ P N . On the other hand, the estimation of the effective channels {ḡ k } K k=1 resulting from covariance shaping requires τ ≥ P and, therefore, the pilot length can potentially be reduced.
• The reference scheme in (B) allows the transmission of up to N data streams per UE, while only one data stream per UE is transmitted when covariance shaping is used. In the following, we demonstrate how the proposed covariance shaping method can effectively outperform spatial multiplexing in scenarios where the UEs exhibit highly overlapping channel covariance matrices and the spatial selectivity of the BS is not sufficient to separate such UEs. We evaluate different scenarios in which the BS serves K closely spaced UEs. We assume that the UEs are equipped with ULAs, whereas both ULAs and uniform planar arrays (UPAs) are considered at the BS. In this context, we define the ULA response at each UE for the angle of impingement θ as where δ = 0.5 represents the ratio between the antenna spacing and the signal wavelength. Likewise, we define the   (37) where M x and M y represent the number of BS antennas along the azimuth and elevation directions, respectively, with M = M x M y . The ULA response at the BS for the angle of impingement φ can be recovered from (37) by setting M x = M , M y = 1, and ψ = 0. Then, the instantaneous channel between the BS and each UE k follows the discrete physical channel model in [38] and is given by where κ is the Ricean factor, U is the number of reflected paths, d k and d k,u are the distances of the LoS path and of the uth reflected path, respectively, θ k and θ k,u are the angles of impingement at the UE of the LoS path and of the uth reflected path, respectively, φ k (resp. ψ k ) and φ k,u (resp. ψ k,u ) are the angles of impingement at the BS of the LoS path and of the uth reflected path along the azimuth (resp. elevation) direction, respectively, α k,u is the random phase delay of the uth reflected path, and β is the pathloss exponent. Unless otherwise stated, we use the simulation parameters listed in Table I. The following results are obtained by means of Monte Carlo simulations with 5 × 10 3 independent channel realizations. Lastly, we point out that, in the considered simulation scenarios, Algorithms 1 and 2 are always observed to converge in very few iterations, i.e., less than 10.

A. Two-UE case
Following the discussion in Section IV-A, we first study the simple case of K = 2 UEs in order to highlight the key features of the proposed covariance shaping framework. In this setting, we assume P = 1, i.e., the UEs are assigned the  Fig. 1(a), with ULA at the BS.
same pilot, and we remark that there is no pilot contamination among the antennas of the same UE for the reference scheme in (B). Assuming that the BS is equipped with a ULA, we consider the LoS and NLoS scenarios illustrated in Fig. 1(a) and in Fig. 2(a), respectively, where the corresponding covariance shaping vectors are shown in Fig. 1(b) and in Fig. 2(b), respectively. Observe that, in both cases, the inter-UE distance is d = 4 m: this is quite small compared with the distance of about 50 m between the UEs and the BS, and results in highly overlapping channel covariance matrices. Let us define the normalized MSE (NMSE) of the channel estimation for UE k as for the effective channels {ĝ k } K k=1 resulting from covariance shaping and as for the channel matrices {Ĥ k } K k=1 . Considering the NLoS scenario in Fig. 1(a), Fig. 3 plots the NMSE of the channel estimation versus the transmit power at the UEs. 6 Thanks to the improved statistical orthogonality between the UEs, covariance shaping considerably increases the channel estimation accuracy with respect to the case where the original MIMO channels are estimated in the presence of pilot contamination among the antennas of different UEs.
Considering again the NLoS scenario in Fig. 1(a), Fig. 4 plots the average sum rate versus the transmit power at the BS ρ BS for different numbers of BS antennas M . For M ≤ 128, it is straightforward to see that the signal subspace separation enforced by covariance shaping during both the uplink pilot-aided channel estimation phase and the downlink data transmission phase has a highly beneficial effect on the system 6 Note that the transmit signal-to-noise ratio (SNR) is given by simply dividing the transmit power at the BS by the noise variance at the BS, which is given in Table I.  performance. However, covariance shaping is outperformed by the reference scheme in (B) for M = 256. In this setting, the enhanced spatial selectivity of the BS allows to effectively separate the UEs and promotes the transmission of multiple data streams per UE. Focusing now on the LoS scenario in Fig. 2(a), Fig. 5

B. Multi-UE Case
Following the discussion in Section IV-B, we now examine the general case of K ≥ 2. We begin by considering the NLoS scenario in Fig. 6(a) with K = 4. Here, the BS is equipped with a ULA and we vary the inter-UE distance d in order to quantify its impact on the degree of statistical orthogonality among the UEs. In this setting, P = 2 orthogonal pilots  are assigned such that two adjacent UEs always utilize orthogonal pilots. The covariance shaping vectors obtained with Algorithm 2 are shown in Fig. 6(b): these tend to focus their power along reflected paths that are as orthogonal as possible to each other while also carrying sufficient channel power (cf. Fig. 1(b)). Fig. 7 plots the average sum rate versus the inter-UE distance d, with M = 128 and ρ BS = 30 dBm. Note that, under the considered pilot assignment, pilot contamination occurs only between UEs at distance 2d. In this context, we also study the case where the UEs are scheduled into two separate groups of non-adjacent UEs without pilot contamination (scheduling), as opposed to the case where all the UEs are the served simultaneously with pilot contamination (no scheduling). As in Fig. 5, the reference scheme in (B) outperforms covariance shaping only for very large inter-UE distances (i.e., d ≥ 13 m). In this regard, the scheduling approach is shown to further deteriorate the performance of spatial multiplexing due to the pre-log factor of 1/2 in the sum rate despite avoiding the pilot contamination. We conclude this section by considering the NLoS scenario in Fig. 8 with K = 8. Here, the BS is equipped with a UPA and is placed at a height of 20 m with respect to the UEs and the scatterers. Moreover, the UEs are separated into two groups with inter-group distance D = 20 m and the inter-UE distance is d = 2 m, which results in highly overlapping channel covariance matrices within the same group. In this setting, P = 4 orthogonal pilots are assigned such that two adjacent UEs always utilize orthogonal pilots and covariance shaping is applied independently to each group. Fig. 9 plots the average sum rate versus the transmit power at the BS ρ BS for different numbers of BS antennas M . Note that M = 64 corresponds to M x = M y = 8, whereas M = 128 corresponds to M x = 16 and M y = 8. As in Fig. 7, we also study the case where the UEs in each group are scheduled into two separate subgroups of non-adjacent UEs without pilot contamination among groups (scheduling). The signal subspace separation enforced by covariance shaping during both the uplink pilotaided channel estimation phase and the downlink data transmission phase becomes even more beneficial when the BS is equipped with a UPA with limited spatial selectivity in the azimuth direction. Again, the scheduling approach is shown to further deteriorate the performance of spatial multiplexing due to the pre-log factor of 1/2 in the sum rate despite avoiding the pilot contamination within the same group.

VI. CONCLUSIONS
In this paper, we introduce the novel concept of MIMO covariance shaping as a means to achieve statistical orthogonality among interfering UEs. It consists in preemptively applying a statistical beamforming at each UE during both the uplink pilot-aided channel estimation phase and the downlink data transmission phase, aiming at enforcing a separation of the signal subspaces of the UEs that would be otherwise highly overlapping. The proposed MIMO covariance shaping framework exploits the realistic non-Kronecker structure of massive MIMO channels, which allows to suitably alter the channel statistics perceived at the BS by designing the transceiver at the UE-side. To compute the covariance shaping strategies, we present a low-complexity block coordinate descent algorithm that minimizes the inter-UE interference (as a metric to measure the spatial correlation), which is proved to converge to a limit point of the original nonconvex problem or to a stationary point in the two-UE case. We provide numerical results characterizing several scenarios where MIMO covariance shaping outperforms a reference schemes employing the multiple antennas at the UE for spatial multiplexing. Specifically, this occurs when the spatial selectivity of the BS is not sufficient to separate UEs exhibiting high spatial correlation, (e.g., when the BS is equipped with a UPA with limited spatial selectivity in the azimuth direction) and the channel estimation is limited by strong pilot contamination. Future work will consider MIMO covariance shaping in combination with multi-stream transmission and will analyze the trade-off between statistical orthogonality and spatial multiplexing.
APPENDIX A DERIVATIONS OF THE EFFECTIVE SINR IN (16) In this Appendix, we derive the expectation terms in (15) to obtain the expression of the effective SINR in (16). First of all, given the multi-UE effective channel matrixH = [ḡ T 1 , . . . ,ḡ T K ] T , where the effective channels {ḡ k } K k=1 are independent andḡ k ∼ CN (0,Φ k ), it is easy to verify that Furthermore, recalling the definition ofĝ k in (9), we have