{'fraction_non_alphanumeric': 0.05167132174069736, 'fraction_numerical': 0.009820704567979097, 'mean_word_length': 3.9494983277591973, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hierarchical Interference Mitigation for Massive MIMO Cellular Networks

We propose a hierarchical interference mitigation scheme for massive MIMO cellular networks. The MIMO precoder at each base station (BS) is partitioned into an inner precoder and an outer precoder. The inner precoder controls the intra-cell interference and is adaptive to local channel state information (CSI) at each BS (CSIT). The outer precoder controls the inter-cell interference and is adaptive to channel statistics. Such hierarchical precoding structure reduces the number of pilot symbols required for CSI estimation in massive MIMO downlink and is robust to the backhaul latency. We study joint optimization of the outer precoders, the user selection, and the power allocation to maximize a general concave utility which has no closed-form expression. We first apply random matrix theory to obtain an approximated problem with closed-form objective. We show that the solution of the approximated problem is asymptotically optimal with respect to the original problem as the number of antennas per BS grows large. Then using the hidden convexity of the problem, we propose an iterative algorithm to find the optimal solution for the approximated problem. We also obtain a low complexity algorithm with provable convergence. Simulations show that the proposed design has significant gain over various state-of-the-art baselines.


I. INTRODUCTION
Massive MIMO is regarded as a promising technology in future wireless networks due to its high spectrum and energy efficiency [1]. The large spatial degree of freedom (DoF) of massive MIMO systems can contribute to (i) spatial multiplexing gains for intra-cell users per BS (MU-MIMO) as well as (ii) inter-cell interference mitigation between the BSs via linear precoders at the BSs. In [2], zero-forcing (ZF) and regularized zero-forcing (RZF) have been proposed for spatial multiplexing of data streams to intra-cell users. More complicated linear precoding schemes based on duality [3] or semidefinite relaxing (SDR) [4] have also been proposed to achieve a better performance. On the other hand, the intercell interference mitigation between BSs is more complicated. One commonly adopted approach to mitigate the inter-cell interference is the coordinated MIMO strategy [5], which performs joint precoding among the BSs using the global realtime CSIs shared among the BSs. Alternatively, cooperative MIMO techniques can also be exploit to mitigate inter-cell interference by sharing both real-time CSI and payload data among the concerned BSs [6].
However, these conventional spatial multiplexing and interference mitigation techniques cannot be applied directly to massive MIMO cellular networks due to the following reasons. First, the MU-MIMO precoding requires real-time local CSIT at the BS. However, the amount of pilot symbols for channel estimation is limited by the coherence time of the channel and it is practically infeasible to obtain good CSI quality when each BS is equipped with a massive MIMO array. Second, the existing inter-cell interference mitigation methods such as cooperative and coordinated MIMO require real-time global CSIT, which is difficult to achieve in practice due to the backhaul latency 1 . Hence, the performance of these schemes is very sensitive to CSIT errors due to outdatedness.
In this paper, we address the above issues by proposing a hierarchical interference mitigation scheme for massive MIMO cellular networks. In the proposed scheme, the MIMO precoder at each BS is partitioned into an inner precoder and an outer precoder as illustrated in Fig. 3. The inner precoder is used to support MU-MIMO (control intra-cell interference and capture the spatial multiplexing gain) at each BS and it is adaptive to real-time local CSIT. The outer precoder can leverage on the remaining spatial DoF to mitigate the inter-cell interference by restricting the transmitted signal at each BS into a subspace and is adaptive to long-term channel statistics 2 . Such hierarchical precoding structure simultaneously resolves both the aforementioned practical challenges. For instance, the issue of insufficient pilot symbols for real-time local CSI estimation is resolved because the BS only needs to estimate the CSI within the subspace determined by the outer precoder, which is of a much smaller dimension than the number of antennas. Furthermore, the outer precoder is adaptive to the long-term channel statistics, which is insensitive to backhaul latency. As a result, the proposed hierarchical precoding framework exploits the spatial DoF to simultaneously achieve spatial multiplexing per BS and inter-cell interference mitigation without expensive backhaul signaling requirement. We consider joint optimization of the outer precoders, the user selection, and the power allocation to maximize a general concave utility function of the average data rates of users. The following first-order challenges need to be addressed.  user set 3 . On the other hand, the optimization of outer precoder also depends on user selection because the outer precoder only needs to suppress the interference to the selected users in other BSs. • Combinatorial Optimization Problem: The user selection problem with hierarchical precoding in the massive MIMO cellular networks is combinatorial with exponential complexity w.r.t. the total number of users. To address the above challenges, we first apply the random matrix theory to obtain an approximated problem with closed-form objective. Then using the hidden convexity of the problem, we propose an iterative algorithm to find the optimal solution for the approximated problem. We also obtain a low complexity algorithm with provable convergence. Finally, we illustrate with simulation that the proposed design achieves significant performance gain compared with various state-ofthe-art baselines under various signaling backhaul latency.
Notations: The superscripts (·) T and (·) † denote transpose and Hermitian respectively. For a set S, |S| denotes the cardinality of S. The operator diag (a) represents a diagonal matrix whose diagonal elements are the elements of vector a. The notation U M ×N denote the set of all M × N semiunitary matrices. Let 1 (·) denote the indication function such that 1 (E) = 1 if the event E is true and 1 (E) = 0 otherwise. span (A) represents the subspace spanned by the columns of a matrix A and orth (A) represents a set of orthogonal basis of span (A). A is the spectral radius of A.

II. SYSTEM MODEL A. Massive MIMO Cellular Network
Consider the downlink of a massive MIMO cellular network with N BSs and K single-antenna users as illustrated in Fig. 1 for N = 2 and K = 5. Each BS has M antennas with M much larger than the number of the associated users. Denote h k,n ∈ C M as the channel vector between BS n and user k. The channel fading process is modeled as h k,n (t) = √ M Θ 1/2 k,n (t) z k,n (t) , ∀k, n, where z k,n (t) ∈ C M has i.i.d. complex entries of zero mean and variance 1/M ; and Θ k,n (t) ∈ C M ×M is the spatial correlation matrix between BS n and user k. The random process z k,n (t) is quasi-static within a time slot but i.i.d. w.r.t. time slots, user and BS indices (t, k, n). The spatial correlation process Θ k,n (t) is assumed to be a slow ergodic process (i.e., Θ k,n (t) remains constant for a large number of time slots) according to a general distribution. As such, the CSI is divided into instantaneous CSI H = {h k,n } and global statistical information Θ {Θ k,n } (spatial correlation matrices). Due to local scattering [7], the spatial correlation matrices of different users in cell n is usually different. However, if the coverage area of a BS is partitioned into N c small sub-areas, it is reasonable to assume that any two users collocated in the same sub-area have almost the same spatial correlation matrices. This motivates us to consider the following locally-clustered spatial channel model. Assumption 1 (Locally-clustered Spatial Channel). The spatial correlation matrices {Θ k,n , ∀k} associated with BS n belongs to a finite set Ψ n with the size |Ψ n | = N c . Furthermore, due to the local spatial scattering [7], we have Rank (Θ k,n ) < M, ∀k, n.
Assumption 1 is realistic because in practice, there are only limited number of significant eigenvalues in a MIMO channel (especially for large M ). The massive MIMO cellular network can be represented by a topology graph as define below.
Definition 1 (Network Topology Graph). For given spatial correlation matrices Θ, define the topology graph of the massive MIMO cellular network as a bipartite graph G T (Θ) = {B, U, E}, where B denotes the set of all BS nodes, U denotes the set of all user nodes, and E is the set of all edges between the BSs and users. For each BS node n, let U n denote the set of associated users and U n = {k : k / ∈ U n , (k, n) ∈ E} denote the set of neighbor users. For each user node k, let b k denote the index of its serving BS and B k = {n : n = b k , (k, n) ∈ E} denote the set of neighbor BSs.
Define E h k,n 2 = Tr (Θ k,n ) as the path gain between BS n and user k. An edge between a user node and a BS node in the network topology graph indicates there is strong path gain between these two nodes. This is stated formally below.
Definition 2 (Edge Set). For given Θ, there is an edge (k, n) ∈ E between BS node n ∈ B and user node k ∈ U in the network topology graph G T (Θ) = {B, U, E} if Tr (Θ k,b k ) < θTr (Θ k,n ), for some threshold θ > 1.
Remark 1. In practical wireless networks, the data rate of each user is limited by the available modulation and coding schemes (MCS) (e.g., the highest MCS in LTE is 64QAM, no coding [8]). If the path gain between a user and a BS is sufficiently small compared to the direct link path gain (θ times smaller than the direct link path gain), the interference from this BS will have negligible effect on the data rate of this user. Simulations show that the performance of the proposed scheme is not sensitive to the choice of θ for a wide range of θ from 5dB to 20dB.
An example of massive MIMO cellular network and the corresponding topology graph is illustrated in Fig. 1. For BS 1, the set of associated users is U 1 = {1, 2}, and the set of neighbor users is U 1 = {3, 4}. For user 2, the index of the serving BS is b 2 = 1 and the set of neighbor BSs is B 2 = {2}. For user 3, the index of the serving BS is b 3 = 2 and the set of neighbor BSs is B 3 = {1}.
At each time slot, linear precoding is employed at BS n to support simultaneous downlink transmissions to a set of scheduled users denoted by S n . Let S = ∪ N n=1 S n denote the set of all the selected users and S n = U n ∩ S denote the set of selected users who are neighbors of BS n. Note that a user k ∈ S n can be potentially interfered by BS n because there is a cross link (edge) between BS n and a user k ∈ S n . For example, consider the massive MIMO cellular network in Since user 3 has a cross link with BS 1 as illustrated in Fig. 1-(a), it can be potentially interfered by BS 1. Using the above notations, the received signal for a user k can be expressed as: where s k ∼ CN (0, 1) is the data symbol, p k is the power allocation and v k is the precoding vector of user k; S b k is the set of selected users at BS b k ; s n = [s l ] l∈Sn ∈ C |Sn| is the data symbol vector at BS n; P n = diag (p n ) and p n = [p l ] l∈Sn ∈ R |Sn| + is the power allocation vector at BS n; V n = [v l ] l∈Sn ∈ C M ×|Sn| is the precoding matrix at BS n; and n k ∼ CN (0, 1) is the AWGN noise.

B. Hierarchical Interference Mitigation
Conventional interference mitigation techniques for small scale MIMO cellular networks such as MU-MIMO precoding, coordinated MIMO [5], or cooperative MIMO [6], cannot be applied directly to massive MIMO cellular networks due to two practical challenges, namely, the insufficient pilot symbols for CSI estimation and the backhaul latency. To resolve these practical challenges, we propose a novel hierarchical interference mitigation control, which can fully utilize the large number of antennas to simultaneously mitigate the intercell interference as well as realize the spatial multiplexing gain per BS. Specifically, the interference mitigation strategy is partitioned into long-term and short-term controls. The shortterm control is responsible to capture the spatial multiplexing gain among the intra-cell users at each BS based on the local CSIT only. On the other hand, the long-term control is responsible to mitigate the inter-cell interference based on the global statistical information. They are elaborated as follows.
1) Hierarchical Precoding for Intra-cell and Inter-cell Interference Mitigation: We first use a simple example to illustrate the idea of hierarchical precoding. Example 1. Consider the massive MIMO cellular network in Fig. 2-(a). Each BS has 32 antennas. Consider the one ring model in [9] for transmit antenna correlation, where a user is surrounded by a ring of scatterers such that the support of its Angle-of-Departure (AOD) distribution are restricted to a certain region as illustrated in Fig. 2-(a). Assume that user 1 and user 2 share the same scattering ring and h 1,1 , h 2,1 are restricted in the same 4-dimensional subspace F (i.e., span (Θ 1,1 ) = span (Θ 2,1 ) = F) as illustrated in Fig. 2-(b). Moreover, due to the local scattering configuration as illustrated in Fig. 2-(a), the cross link channel vector h 3,1 between BS 1 and user 3 is restricted in a 2-dimensional subspace F c ⊂ F (i.e., span (Θ 3,1 ) = F c ⊂ F). We consider a hierarchical precoder structure for BS 1: where F 1 ∈ U 32×2 with span (F 1 ) ⊆ F is the outer precoder adaptive to the spatial correlation matrices Θ only, and G 1 ∈ C 2×2 is the inner precoder adaptive to the local CSI (the effective channel We have the following observations from Example 1. Role of spatial correlation matrices Θ: The knowledge of the spatial correlation matrices Θ can be exploited to design the outer precoder F 1 to eliminate the inter-cell interference to user 3. Specifically, this can be achieved by setting F 1 to be a set of orthogonal basis of a 2-dimensional subspace F 1 , where F 1 ⊂ F and F 1 ⊥ F c as illustrated in Fig. 2-(b).
Role of local CSI H 1 : The knowledge of local real-time instantaneous CSI H 1 can be exploited to design the inner precoder G 1 to realize the spatial multiplexing gain at BS 1.
For general massive MIMO cellular networks, we propose a hierarchical precoder V n = F n G n for each BS n as illustrated in Fig. 3. The outer precoder F n ∈ U M ×Mn with M n < M (we let F n = 0 if M n = 0) is used to eliminate the inter-cell interference and is adaptive to the global statistical information Θ. The inner precoder G n ∈ C Mn×|Sn| is used to realize the spatial multiplexing gain at each BS and is adaptive to the local real-time CSI H Sn F † n h k,n † k∈Sn ∈ C |Sn|×Mn . Define F = {F 1 , ..., F N } as the set of outer precoders for all BSs. By properly choosing F (equation (2)), one can eliminate the inter-cell interference as illustrated in Example 1. Remark 2. Physically, the rank M n of F n means the number of data streams for spatial multiplexing at BS n. Due to limited spatial scattering [7], a BS with say M = 100 antennas does not mean it can support spatial multiplexing of 100 data streams. In practice, there are just a few (say 10) significant eigenchannels despite having 100 antennas, and having M n = 10 spatially multiplexed data streams already capture most of the spatial multiplexing advantage. The remaining spatial DoFs can be used for inter-cell interference mitigation.
For a given outer precoder F, we consider regularized zeroforcing (RZF) inner precoder with a parameter ν. The RZF precoder is easy to implement and is asymptotically optimal for M, |S n | → ∞ [10]. Moreover, we can apply the technique of deterministic equivalent (DE) for RZF precoding in [11] to facilitate the algorithm design. For convenience, define the composite channel from BS n to any subset of users U S ⊆ U as H U S = [h l,n ] † l∈U S ∈ C |U S |×M . If the inter-cell interference is completely eliminated by the outer precoders F, the RZF inner precoder is given by 2) Statistical User Scheduling and Power Allocation: In [12], it has been observed that as the number of antennas M grows large, the role of multi-user diversity gain (by selecting users based on instantaneous CSIT) becomes less and less effective because of "channel hardening". Moreover, the benefits of short timescale power allocation (i.e., the power allocation is adaptive to instantaneous CSIT) becomes asymptotically negligible as K, M → ∞ because the data rate of each user converges almost surely to a deterministic function of the power allocation vector p = p T 1 , ..., p T N T as will be shown in Lemma 1. As such, the user selection set and power allocation is assumed to be adaptive to the global statistical information Θ only. Specifically, the user selection and the outer precoder are chosen to satisfy the zero inter-cell interference constraint: and the power allocation has to satisfy a per-BS power constraint that will be elaborated later. Note that the inter-cell interference from BS n to a user k ∈ S n can be expressed as Hence, the constraint in (2) can be used to eliminate the inter-cell interference to all users in the system. In (2), we consider ZF criteria for inter-cell interference mitigation due to its simplicity and asymptotic optimality at high SNR. Similar ZF criteria has also been used in [9] to design pre-beamforming matrix based on spatial correlation matrices in single cell systems.
Remark 3 (Implementation Considerations). In the proposed hierarchical interference mitigation, the long term controls such as outer precoding, statistical user selection and power allocation are implemented at a central node based on the global statistical information about the channel (Θ), while the short term control (inner precoding) is implemented locally at each BS based on the local real-time instantaneous CSI knowledge ( H Sn ). The proposed hierarchical precoding solution has several unique benefits regarding implementation. (a) Robust to CSI signaling latency in backhaul, (b) Resolve the issues of insufficient pilot and feedback overhead for real-time CSIT estimation in massive MIMO systems. For instance, the central node requires the global statistical information (spatial correlation matrices) Θ to compute the long term controls {F, S, p}. The spatial correlation matrices can be estimated via downlink training using some standard covariance matrix estimation technique [13] at the users and then fed back to the BSs. Since Θ changes at a much slower time scale w.r.t. the time slot rate, such a design requires substantially less signaling overhead compared to the coordinated MIMO and is more robust w.r.t. to the backhaul latency. On the other hand, BS n only needs to know the local real-time instantaneous CSI H Sn ∈ C |Sn|×Mn for the inner precoder G n . This can be obtained via downlink channel estimation and channel feedback using conventional CSI signaling mechanisms in modern wireless systems such as LTE [14]. Since M n is substantially smaller than M , the issue of the huge downlink pilot and CSI feedback signaling overhead in massive MIMO is also alleviated by the hierarchical precoding design.

III. OPTIMIZATION FORMULATION FOR HIERARCHICAL INTERFERENCE MITIGATION
We consider joint optimization of the outer precoders F, the user selection S, and the power allocation p = p  (2), the instantaneous data rate (treating interference as noise) of user k is where the precoders [v l ] l∈Sn = F n G n , ∀n with the inner precoder G n given by (1). The transmit power of BS n is Note that there may not always be enough spatial DoFs to eliminate the inter-cell interference to all the users. Hence, for a fixed composite control variable Γ, it is possible that only part of the users can be scheduled for transmission. For fairness considerations, we consider randomized control policy which realizes time-sharing between several composite control variables as defined below.
Definition 3 (Randomized Control Policy). A randomized control policy Ω = {Ξ, q} consists of a set of composite control variables Ξ Γ 1 , ..., Γ |Ξ| with |Ξ| ≤ K and a probability vector q At any time slot, the composite control variable Γ j is used with probability q j , i.e., the outer precoders, the user selection set and the power allocation are respectively given by F (j), S (j) and p (j) with probability q j . Moreover, define the set of feasible control policies under per-BS power constraint P c as For given control policy Ω = {Ξ, q} and spatial correlation matrices Θ, the conditional average data rate of user k is: The network performance is characterized by a utility function U (r (Ω|Θ)), where r (Ω|Θ) = [r 1 (Ω|Θ) , ..., r K (Ω|Θ)] T is the conditional average rate vector. We make the following assumptions on U (r) (r is a simplified notation for r (Ω|Θ)).
Assumption 2 (Assumptions on Utility). The utility function can be expressed as U (r) K k=1 w k u (r k ), where w k ≥ 0 is the weight for user k, u (r) is assumed to be a twice differentiable, concave and increasing function for all r ≥ 0. Moreover, u (r) is L-Lipschitz, i.e., The above utility function captures a lot of interesting cases: • Alpha-fair [15]: Alpha-fair can be used to compromise between the fairness to users and the utilization of resources. The utility function is 4 where > 0 is a small number. • Proportional Fair (PFS) [16]: This is a special case of alpha-fair when α = 1. For a given topology graph G T (Θ) = {B, U, E} and per-BS power constraint P c , the problem of interference mitigation via hierarchical precoding can be formulated as 5 : Note that the conditional average rate r k (Ω|Θ) in the utility function and the conditional average power E [ P n (Γ)| Θ] in the constraint function of P (G T (Θ)) do not have closed form expressions. To make the problem tractable, we need to address the following challenge.
Find an approximated problem P E (G T (Θ)) with closed form objective and constraints such that the solution of P E (G T (Θ)) is asymptotically optimal w.r.t. P (G T (Θ)) as M grows large.
We resort to random matrix theory to solve the above challenge. Specifically, we first derive deterministic equivalents (DEs) [11] for the conditional average rate and power. Then we obtain an approximated problem P E (G T (Θ)) by replacing the conditional average rate and power with their DE approximations. Finally, we show that the solution of P E (G T (Θ)) is an O (ν)-optimal solution of P (G T (Θ)). if it satisfies the zero inter-cell interference constraint F † n (j) k∈Sn(j) U n∩S (j) Θ k,n = 0, ∀n, j and the following relaxed per-BS power constraint where U * is the optimal objective value of P (G T (Θ)).
Assumption 3-1) is satisfied by many MIMO channel models such as the angular domain MIMO channel model in [7] and it is a standard assumption in the literatures, see e.g., [11], [17]. Under Assumption 3-1), Assumption 3-2) holds true if lim sup M →∞ |{Θ k,n : k ∈ U n }| < ∞, that is, if {Θ k,n : k ∈ U n } belongs to a finite family [11]. According to the locally-clustered channel model in Assumption 1, we have|{Θ k,n : k ∈ U n }| ≤ N c and thus Assumption 3-2) holds true. Assumption 3-3) is to ensure that the utility function is bounded as K → ∞.
are the deterministic equivalent (DE) of user rate and BS transmit power, and ξ i , ∀i ∈ S n form the unique solution of withΘ i,n = F n F † n Θ i,n F n F † n , ∀i ∈ S n . Please refer to Appendix A for the proof. Remark 4. Note that the above DEs are established on the "conditional distribution of the channel" (conditioned on the statistics Θ). Given a realization of Θ (the statistics), the control actions Γ = {F, S, p} are all fixed (because they are adaptive to Θ only). As such, the conditional measure of H (conditioned on the given Θ) will exhibit "random matrix theory" behavior and the DE convergence in Lemma 1 can be proved using standard techniques in [11]. On the other hand, if S were adaptive to the instantaneous CSI H (short-term user selection), then conditioned on Θ, S would be random and hence the DE approximation would fail (due to the random S or extreme value effect of the user selection which changes the underlying conditional distribution of the channels H). Similar conclusion has also been made in [9] that the DE of the data rate in massive MIMO system is valid as long as the user selection is independent of the instantaneous CSI H.
Based on Lemma 1, we have the following result.
Let Ω denote the optimal solution of Given Assumption 3 and for sufficiently small ν > 0, Ω is an O (ν)-optimal solution of P (G T (Θ)) as M → ∞.
Please refer to Appendix B for the proof. By Theorem 1, the solution of P (G T (Θ)) can be approximated by the solution of P E , and the approximation is O (ν)-optimal as M → ∞.

IV. SOLUTION OF P E (G T (Θ))
In the rest of the paper, we focus on solving P E (G T (Θ)) for fixed Θ. We will use G T , r • k (Γ) and P • n (Γ) as simplified notations for G T (Θ), r • k (Γ|Θ) and P • n (Γ|Θ) when there is no ambiguity. Clearly, the utility function U E (Ω) is not a convex function on Ω and thus P E (G T ) is a non-convex optimization problem. Moreover, the optimization variables in P E (G T ) involve a set of composite control variables Ξ with undetermined size and the associated probabilities q with undetermined dimension. It is in general very difficult to find the global optimal solution for such a non-convex problem. In this section, we are going to address the following challenge.
Challenge 2 (Design a Global Convergent Algorithm for P E (G T )). Exploit the specific structure of problem P E (G T ) to design an iterative algorithm that converges to the global optimal solution of P E (G T ).
We first study the optimality condition of P E (G T ). Then we propose an iterative algorithm to solve Challenge 2.

A. Global Optimality Condition of P E (G T )
It is difficult to find a simple characterization for the necessary and sufficient global optimality condition of a general non-convex problem. However, problem P E (G T ) is not an arbitrary non-convex problem but has some specific hidden convexity structure, which can be exploited to derive the global optimality condition for P E (G T ) as shown below.
We first study the hidden convexity of P E (G T ). Define the (deterministic equivalent of) average rate region as: where r • (Ω) = [r • 1 (Ω) , ..., r • K (Ω)] T with r • k (Ω) = |Ξ| j=1 q j r • k (Γ j ). Then we have the following Lemma. Lemma 2 (Convexity of R). R = Conv R F , where Conv (·) denotes the convex hull operation and R F r Please refer to Appendix C for the proof.
The following lemma shows that problem P E (G T ) is equivalent to a convex problem: Lemma 3 (Equivalence between P E (G T ) and (11)). If Ω is the global optimal solution of P E (G T ), then r • (Ω ) is the optimal solution of problem (11); on the other hand, if r • is the optimal solution of problem (11), then any Ω satisfying r • (Ω ) = r • is also the global optimal solution of P E (G T ).
Please refer to Appendix C for the proof. This hidden convexity of P E (G T ) (i.e., the equivalence between P E (G T ) and (11)) is the key to derive the global optimality condition of P E (G T ). Note that although problem (11) is convex, the solution is still non-trivial because there is no simple characterization for its feasible set R.
To derive the global optimality condition of P E (G T ), we also need the first order optimality condition of problem (11) as summarized in the following lemma.
Finally, from Lemma 3 and Lemma 4, we can obtain the necessary and sufficient global optimality condition for problem P E (G T ) as follows.

B. Global Optimal Solution of P E (G T )
Just as we can obtain the optimal solution of a convex problem by solving its KKT conditions, we can also obtain the global optimal solution of P E (G T ) by solving the global optimality condition in Theorem 2. Specifically, for any given spatial correlation matrices Θ, we propose Algorithm E to achieve the global optimality condition of P E (G T (Θ)) by iteratively updates the optimization variables Ξ, q and the weight vector µ in Theorem 2. Algorithm E (Top level algorithm for solving P E (G T (Θ))): Initialization: Set i = 0 and let µ (0) = [w k ] k=1,...,K . Call Procedure W with input µ (0) to obtain a composite control variable Γ µ (0) and let Ξ (0) = Γ µ (0) .
Step 1 (Update probability vector q): Let Step 2 (Update composite control variable set Ξ): Let (14) Call Procedure W with input µ (i+1) to obtain a new composite control variable Γ µ (i+1) . Update Ξ as Step 3: If i > 0 and UE Ω (i) − UE Ω (i−1) ≤ ε, where ε > 0 is a small number, terminate the algorithm. Otherwise, let i = i + 1 and return to Step 1. Fig. 4 summarizes the inter-relationship between the components of Algorithm E. Algorithm E contains two procedures (subroutines) which will be elaborated below. Remark 5. Algorithm E can be interpreted as the Frank-Wolfe Algorithm (also known as the conditional gradient algorithm) with exact line search [18] applied on the equivalent convex problem in (11). Compared to the conventional Frank-Wolfe Algorithm, the main difference is that the optimization variable in problem P E (G T ) is Ω instead of r • in (11), and the optimization w.r.t. Ω is non-convex. Nonetheless, we can exploit the hidden convexity (global optimality condition) in Theorem 2 to establish the global convergence of Algorithm E as will be shown in Theorem 4.
1) Procedure Q (Optimization of q for fixed Ξ): For given input Ξ, Procedure Q essentially solves the optimal probability vector q for P E (G T ) with fixed Ξ, i.e., Procedure Q with input Ξ is a standard convex optimization procedure to solve the following optimization problem: s.t. q j ∈ [0, 1] , ∀j and |Ξ| j=1 q j = 1, where Γ j is the j-th composite control variable in Ξ. Hence, Procedure Q can be efficiently implemented by existing convex optimization methods/software. As such, the pseudo code of Procedure Q is omitted here for conciseness.
2) Procedure W (Finding a new composite control variable for given µ): The pseudo code of Procedure W is summarized in Table I. In Line 2, ξ i , ∀i ∈ S n is the unique solution of (9) withΘ i,n =Θ i,n S n where U S n = orth k∈Sn Θ k,n . In Line 4, R (S) is the (deterministic equivalent of) weighted sum-rate for given user selection S. In Line 7, S n = U n ∩ S . For convenience, ξ k is referred to as the effective channel gain of user k andΘ i,n S n is called the projected spatial correlation matrix of user i. For conciseness,Θ i,n S n is denoted as Θ i,n when there is no ambiguity. To calculate the weighted sum-rate R (S), we need to obtain the effective channel gains ξ i , ∀i ∈ S associated with S by solving the fixed point equation in (9). The solution of (9) can be obtained using the following fixed point iterations [11] For given input µ, Procedure W essentially finds a composite control variable Γ (µ) which satisfies the global optimality condition in (13) for fixed µ.
Theorem 3 (Characterization of Procedure W ). For given input µ, the output Γ (µ) of Procedure W satisfies Please refer to Appendix D for the proof.
3) Convergence and Performance of Algorithm E: The update rule in Algorithm E is designed according to the global optimality condition in Theorem 2. As a result, it can be shown that Algorithm E converges to the global optimal solution of P E (G T ) using the global optimality condition in Theorem 2 and the property of Algorithm E in the following Lemma.

Lemma 5 (Property of Algorithm E).
Let Ω (i) be the control policy in the i-th iteration of Algorithm E. We have where µ (i+1) = ∇U r • Ω (i) is given in (14).
Please refer to Appendix E for the proof. Using Theorem 2 and Lemma 5, we obtain the following global convergence result.
Theorem 4 (Global Optimality of Algorithm E). Algorithm E monotonically increases the utility U E Ω (i) and lim i→∞ U E Ω (i) → U , where U is the global optimal value of P E (G T ).
Please refer to Appendix E for the proof. In step 2 of Algorithm E, we need to call Procedure W , which involves an exhaustive user selection process where R (S) is calculated for all possible user set S (see Line 1 to Line 6 of Procedure W ). The complexity of exhaustive user selection is exponential w.r.t. the number of users K. In the next subsection, we will propose a low complexity solution, named the modified Algorithm E, for P E (G T ) by replacing the exhaustive user selection process with a statistical greedy user selection process.

C. Low Complexity Solution of P E (G T )
The low complexity solution (modified Algorithm E) is obtained by replacing the exact solution Γ (µ) of (18) in step 2 (and the initialization step) of Algorithm E with an approximate solutionΓ (µ) found by a low complexity procedure named Procedure W. In other words, the modified Algorithm E are the same as Algorithm E except that Procedure W (which involves exhaustive user selection) is replaced by the low complexity counterpart Procedure W (which is based on statistical greedy user selection).
The pseudo code of Procedure W is summarized in Table  II. In Line 3 and 4, the weighted sum-rate R (S) for any given S can be calculated using the same method as described in Procedure W . Clearly, the statistical greedy user selection loop between Line 2 and Line 9 converges to a solutionŜ within K iterations. Fig. 4 summarizes the overall low complexity solution and the inter-relationship between the components of the modified Algorithm E. To justify the modified Algorithm E, we need to address the following challenge.
Challenge 3 (Monotone Convergence of the modified Algorithm E). Prove the monotone convergence of the modified Algorithm E as well as characterize the performance loss of the modified Algorithm E w.r.t. the global optimal solution.
Theorem 5 (Convergence of the Modified Alg. E). The modified Algorithm E monotonically increases the utility U E Ω (i) and lim i→∞ U E Ω (i) →Û . Moreover, the gap ofÛ with the optimal utility U E of P E (G T ) is bounded by where r • k ,μ can be any accumulation point of the iterates r • Ω (i) , µ (i+1) generated by the modified Algorithm E, Γ (μ) is the output of Procedure W with inputμ andΓ (μ) is the output of Procedure W with inputμ.
Please refer to Appendix F for the proof. Theorem 5 states that the performance gap between the modified Algorithm E and (the optimal) Algorithm E is upper bounded by the performance gap (in terms of weighted sum-rate) between Procedure W (statistical user selection) and Procedure W (exhaustive user selection).
Complexity Analysis for the Modified Algorithm E: The computation complexity is evaluated in terms of the number of matrix multiplications, matrix inversions and Gram-Schmidt processes, since these operations dominate the first order of the overall computation complexity. For simplicity, we assume Rank (Θ i,n ) = d, ∀i, n and U n = K, ∀n. Suppose that the fixed point iterations in (17) converges to the desired accuracy in C f iterations. Then the complexity of Procedure W is analyzed as follows. In each iteration, the greedy search to find the S • requires evaluating K − |S| < K weighted sumrates R S . Each R S needs no more than O (C f |S|) M × M matrix multiplications and O (N C f ) M × M matrix inversions. The greedy search is repeated for at most K times. For each n,Θ i,n , ∀i ∈ U n is updated for at most KK times, and each update needs one M × M Gram-Schmidt processes (i.e., Gram-Schmidt processes for a M × M matrix) to calculate orth

V. SIMULATION RESULTS
Consider a cellular network with 19 cells as illustrated in Fig. 5. The inter-site distance is 500m. In each cell, there are 2 uniformly distributed hotspots with a radius of 50m. There are 12 users in one cell, 2/3 of whom are clustered around the hotspots, while the others are uniformly distributed within the cell. Each BS is equipped with M = 48 antennas. The spatial correlation matrices are generated according to Θ k,n = L k,nΘk,n , ∀k, n, where the path gains L k,n 's are generated using the path loss model ("Urban Macro NLOS" model) in [19], and the normalized spatial correlation matricesΘ k,n 's with Tr Θ k,n = M and Rank (Θ k,n ) = 6 are randomly generated. In the simulations, we set the threshold in Definition 2 as θ = 10dB and the parameter for RZF as ν = 10 −2 . We compare the performance of the proposed algorithm with the following two baselines.
Baseline 1 (FFR): Fractional frequency reuse (FFR) [20] is applied to suppress the inter-cell interference. In each cell, ZF beamforming is used to serve the users on each subband.
Baseline 2 (Clustered CoMP): 3 neighbor BSs form a cluster and employ cooperative ZF [21] to simultaneously serve all the users within the cluster.

A. Convergence of the Modified Algorithm E
Consider the PFS utility U (r) = 1 K K k=1 log (r k + ) with = 10 −4 . The per BS transmit power is P c = 10dB. In Fig. 6, we plot the objective value U E (Ω) of P E versus the number of iterations of the modified Algorithm E. It can be seen that the modified Algorithm E quickly converges.

B. Performance Evaluation under PFS Utility
The simulation setup is the same as that in Fig. 6. In Fig. 7, we compare the average cell throughput of different schemes under different backhaul latencies. For baseline 2, the 3 cooperative BSs need to exchange CSI and payload data, and thus there is CSI delay when the backhaul latency is not zero. When there is CSI delay, the outdated CSI is related to the actual CSI by the autoregressive model in [22]. It can be seen that the cell throughput of the proposed scheme is close to the baseline 2 with zero backhaul latency and is much larger than baseline 1. The worst 10% users also benefit from huge throughput gain   Table III compares the computational complexity (CPU time) and signaling overhead of different schemes. The computational complexity and the backhaul signaling overhead of the proposed scheme are similar to FFR, and are much lower than CoMP. The real-time CSI estimation overhead of the proposed scheme is lower than both FFR and CoMP.

C. Performance Evaluation under Sum-rate Utility
Consider the sum-rate utility. In Fig. 8, we plot the average cell throughput K N U (r) of different schemes versus the per BS transmit power P c . It can be seen that the cell throughput of the proposed scheme is close to the baseline 2 with zero backhaul latency and is much larger than baseline 1. When there is a backhaul latency of 10ms, the proposed scheme also has a significant throughput gain over baseline 2. The DE of   Pc. The user speed is 3 km/h. the cell throughput K N U E (Ω) is also plotted for the proposed scheme. It can be seen that the DE is very accurate.

VI. CONCLUSION
We propose a hierarchical interference mitigation scheme for massive MIMO cellular networks. The MIMO precoder is partitioned into inner precoder (for intra-cell interference control) and outer precoder (for inter-cell interference control). We study joint optimization of the outer precoders, the user selection, and the power allocation. The optimization only requires the knowledge of spatial correlation matrices and thus is robust to backhaul latency. We first apply the random matrix theory to obtain an approximated problem which is non-convex. Then using the hidden convexity of the problem, we propose Algorithm E to obtain the global optimal solution and a low complexity version of Algorithm E to find a suboptimal solution. Simulations show that the proposed design achieves significant performance gain over various state-ofthe-art baselines.
For composite control variable Γ satisfying the conditions in Lemma 1, it can be shown that r k (Γ) and P n (Γ) are uniformly integrable [23] w.r.t. M . Together with Lemma 1, it follows that By definition, we have Then it follows from (26) and (27) that Similarly, it can be shown that We expand U r Ω * |Θ − U r Ω |Θ as follows From (25), we have for Ω ∈ Ω * ,Ω . Then it follows from (31) and From (24,29) and the definition of Ω and Ω * , we have Then it follows from (24,30,32,33) that This completes the proof for Theorem 1.

C. Proofs for the Results in Section IV-A
Proof of Lemma 2: Clearly, R ⊆ Conv R F . Hence, we only need to prove that any Pareto boundary point r of Conv R F must lie in R. First, it is easy to see that r can always be expressed as a convex combination of K + 1 points {j : q j > 0} must lie in the supporting hyperplane to Conv R F at the Pareto boundary point r . Otherwise, r cannot be a Pareto boundary point of Conv R F . The above two facts imply that r can be expressed as a convex combination of K ≤ K points in the set {r • (Γ 1 ) , ..., r • (Γ K+1 )}, i.e., r = K j=1 q j r • Γ j , where K j=1 q j = 1, q j ∈ [0, 1] and Γ j ∈ {Γ 1 , ..., Γ K+1 }. Hence, r must lie in R.
Proof of Lemma 3: The first part of Lemma 3 follows directly from the definition of problem (11) and P E (G T ). The second part of Lemma 3 can be proved by contradiction.
Suppose Ω satisfies r • (Ω ) = r • but is not the global optimal solution of P E (G T ). Then there exists a control policy Ω ∈ Λ • (P c ) such that U E (Ω) > U E (Ω ). Then compared to r • , r • (Ω) ∈ R achieves a larger objective value for problem (11), which contradicts with the assumption that r • is the optimal solution of problem (11).
By Lemma 4, r • (Ω ) is the optimal solution of problem (11). Then it follows from Lemma 3 that Ω is the global optimal solution of P E (G T ).
On the other hand, suppose Ω is the optimal solution of P E (G T ). By Lemma 3, r • (Ω ) is the optimal solution of (11). Then by Lemma 4, Ω satisfies (36), from which it can be shown that Ω satisfies the optimality condition in (13).

D. Proof of Theorem 3
It can be seen that the optimal solution of the following WSRM problem satisfies (18) Hence, we only need to prove that the output Γ (µ) of Procedure W is the optimal solution of P W (G T , µ). First, we show that P W (G T , µ) is equivalent to a joint user selection and power allocation problem. Then Γ = {F , S , p } is an optimal solution of P W (G T , µ), where F = {F 1 , ..., F N } with F n = orth I M − U S n U † S n k∈S n Θ k,n ; and S n = U n ∩ S .
Proof: Lemma 7 can be proved by contradiction. First, it is easy to see that Γ = {F , S , p } is a feasible solution of P W (G T , µ), i.e., Γ ∈ Ξ F• (P c ). Suppose that Γ is not an optimal solution of P W (G T , µ). Then there exists Γ = {F, S, p} ∈ Ξ F• (P c ) such that K k=1 µ k r • k (Γ) > K k=1 µ k r • k (Γ ). Since F satisfies the zero inter-cell interference constraint in (2), we must have span F n F † n ⊆ span I M − U S n U † S n , ∀n. Let Γ = F, S, p , where F = F 1 , ..., F N with F n =orth F n F † n k∈Sn Θ k,n . It can be shown that r • k Γ = log (1 + p k ) = r • k (Γ) , ∀k and P • n Γ = P • n (Γ). Let Γ = F , S, p , where F = F 1 , ..., F N with F n = orth I M − U S n U † S n k∈Sn Θ k,n . It is easy to see that Γ satisfies (2) and r • k Γ = log (1 + p k ) = r • k (Γ) , ∀k. Using the fact that span F n F † n ⊆ span I M − U S n U † S n , it can be shown that P • n Γ ≤ P • n Γ ≤ P c , which implies that S, p is a feasible solution of Problem (38). Hence, we have K k=1 µ k r • k (Γ ) ≥ K k=1 µ k r • k Γ = K k=1 µ k r • k (Γ), which contradicts with K k=1 µ k r • k (Γ) > K k=1 µ k r • k (Γ ). This completes the proof.
It can be verified that S , p in Line 6 of Procedure W is the optimal solution of (38). By Lemma 7, the output Γ (µ) of Procedure W is the optimal solution of P W (G T , µ).
Proof of Theorem 4: Using the fact that any Pareto point of a K-dimensional convex polytope in R K + can be expressed as a convex combination of no more than K vertices, it can be shown that there are at most K non-zero elements in q (i) in step 1 of Algorithm E. Hence Ξ (i) ≤ K, ∀i and the solution found by Algorithm E is feasible.

Lemma 8.
Let Ω (i) be the iterates generated by Algorithm E. We have lim i→∞ U E Ω (i) → U * for some U * .
By Assumption 2, u (r) is L-Lipschitz, which implies that U (x) is also L-Lipschitz with the "L constant" given by L ≤ L max k w k . It is well know that the following lemma holds for a L-Lipschitz function.