the same physical dimension, which allows for large-scale spatial multiplexing and highly directional beamforming. For massive multiple-input-multiple-output (MIMO) mmWave systems, the conventional fully-digital (FD) precoding scheme, which requires one radio frequency (RF) chain per antenna element, is not viable due to the high cost and the high power consumption of RF chain components in high frequencies. To address the challenge of this hardware limitation while exploiting the multiplexing gain of MIMO, hybrid precoding architectures have been proposed for mmWave systems and widely investigated in the literature [6]- [19].

The key idea of hybrid precoding is using a linear network of variable phase shifters in the RF domain in addition to the baseband digital precoding. Such a scheme was firstly known as soft antenna selection (SAS) which was proposed in [20] to improve the performance of the traditional antenna selection scheme. It was shown in [20] that the SAS method can even achieve the optimal performance of a fully-digitally precoded single-stream MIMO system when each end is equipped with two or more RF chains. Recently, SAS is re-introduced to mmWave system design under the name of hybrid precoding/beamforming [2], [6], which has received significant interest in recent literature on large-scale antenna array systems.

The pioneering work on hybrid precoding for mmWave systems [6] investigated hybrid precoding methods for a point-topoint mmWave MIMO system. It was shown in [6] that the spectral efficiency maximization problem for mmWave MIMO systems can be approximately solved by addressing a matrix approximation or reconstruction problem, i.e., minimizing the Frobenius norm of the difference between the optimal FD precoder and the hybrid precoder. By exploiting the sparse nature of mmWave channels, the matrix approximation problem is reformulated as a compressive-sensing-like problem which is addressed by using a modified orthogonal matching pursuit (OMP) algorithm. The OMP method in [6] can achieve good performance when hundreds of antennas are used at both ends of transceiver and/or the number of RF chains is strictly greater than the number of data streams. But there is still a significant performance gap between the hybrid precoding method proposed in [6] and the optimal FD precoding method especially when the number of RF chains is equal to the number of data streams. Hence, several works have devoted to reducing this performance gap. The authors of [7] proposed using alternating optimization to approximate the solution of the matrix approximation problem. To deal with the difficulty arising from the unitmodulus constraints, Yu et al. [7] used manifold optimization to solve the analog (or RF) precoder design problem with fixed digital precoder. Other works that are based on matrix approximation can be found in [8]- [10]. Together with [6] and [7], all the above methods don't directly take into consideration the power constraint in the optimization, they inevitably incur some performance loss as compared to the optimal hybrid precoding performance. Differently from the matrix approximation methods [6]- [10], the work [15] proposed another hybrid precoding algorithm that can approximately solve the spectral efficiency maximization problem of mmWave MIMO systems. First, by removing the coupling among the hybrid precoder and decoders (or combiners) using some approximations under large-antenna array, the analog precoder is designed independently. Then, given the analog precoder, the digital precoder, the analog combiner, and the digital combiner are sequentially designed. Simulations confirm that the hybrid precoding method in [15] can achieve a performance that is very close to the FD precoding performance, though it is still a heuristic algorithm.

While the above works considered hybrid precoding for single-user MIMO systems, there are a few works on hybrid precoding for spectral efficiency optimization of multi-user MIMO/MISO systems [11], [12], [15]- [18]. In [11], the authors developed a low-complexity two-stage hybrid precoding algorithm for downlink mmWave single-stream multi-user MIMO systems with limited feedback channel. In the first stage, the analog precoder at the base station (BS) and the analog combiners at the users are jointly designed to maximize the desired signal power of each user while neglecting the resulting interference among users. In the second stage, the BS digital precoder is designed to manage the multi-user interference. The authors of [12] proposed equal-gain-transmission-based analog beamforming combined with block diagonalization (BD) based digital precoding for generic channel setup. By exploiting the knowledge of angle of departure, the work [13] proposed iterative matrix decomposition based BD for mmWave multiuser MIMO systems. In [14], the authors first proposed a low complexity mmWave channel estimation algorithm for multiuser mmWave systems with single RF chain at each user and then designed a hybrid analog precoding and digital zero-forcing precoding scheme based on the channel estimations. In [15], assuming perfect channel state information, the authors proposed an iterative hybrid precoding algorithm for multi-user MISO systems, which iterates between the design of zero-forcing (ZF) precoder and the analog precoder. The works [16], [17] proposed low complexity hybrid precoding schemes for multi-user MISO systems by performing phase-only maximum ratio combining (MRC) or channel matched filter in the analog domain and ZF beamforming in the digital domain based on the effective channel. Extending the hybrid precoding methods in [16] and [17], the work [18] proposed similar low complexity hybrid precoding algorithms for multi-user MIMO systems by designing the analog and digital precoders sequentially, where digital precoders are separated to pre-filters and post-filters. First, analog precoders and decoders are designed based on a per-user channel matching criterion. Then digital pre-filters are applied at both the transmitter side and the receiver side. Finally, optimal digital post-filters are derived based on four linear transmit strategies, including ZF [21], MMSE [22], block diagonalization (BD) [23], regularized BD (RBD) [24], followed by optimal power allocation via the water-filling algorithm. The simulation results in [18] show that the RBD method performs the best among various transmit strategies. Differently from the above works focusing on the downlink systems, the authors of [19] investigated the hybrid beamforming design problems for a mmWave massive multicell MIMO uplink transmission system, tackling both the intracell and inter-cell interference under the hardware constraints.

For tractability, most of the above works typically assume infinite resolution phase shifters used for analog precoders, although it is very expensive to realize accurate phase shifters [25], [26]. Usually, low resolution but less expensive phase shifters are used in practice, resulting in spectral efficiency optimization problems with discrete unit modulus constraints [6], [15], [17]. Since it is impossible to directly solve the discrete optimization problem, the most straightforward way to design analog precoder with finite resolution phase shifters is to address the spectral efficiency problem assuming infinite resolution first and then to quantize analog precoder using a finite set of available phases [17]. However, this approach could be ineffective for the case of very low resolution phase shifters. An alternative way to approximately addressing the discrete spectral efficiency problem was proposed in [15], where the discrete unit modulus constraints are directly taken into consideration in the optimization of analog precoder.

To the best of our knowledge, all the existing hybrid precoding methods are heuristic and there is no work that tries to directly solve the spectral efficiency optimization problems of various MIMO systems under hybrid precoding setups from the perspective of mathematical optimization. While some of the existing hybrid precoding methods can achieve very high spectral efficiency when the number of RF chains is strictly greater than the number of symbols and even the same performance as the fully-digital precoding method when the number of RF chains is no less than twice the number of streams, there is in some cases a significant gap between the achievable rate of the existing methods and the theoretical maximum rate. Hence, an interesting and fundamental question is: if the spectral efficiency optimization problem is directly addressed from a perspective of mathematical optimization, is it possible for the hybrid precoding method to obtain a better system spectral efficiency and how close it is to the fully-digital precoding performance when the minimum number of RF chains are used. To answer this fundamental question, we propose an iterative hybrid precoding algorithm for multi-stream multi-user MIMO systems, aiming to directly solve the spectral efficiency optimization problem. Apparently, the main difficulty of the spectral efficiency optimization problem arises from the unit modulus constraints and the coupling among the analog precoder and the digital precoders. To address these difficulties, we apply penalty dual decomposition (PDD) method [27], [28] to the spectral efficiency optimization problem by first penalizing and dualizing the coupling constraint into the objective, and then iteratively solving the augmented Lagrangian problem using block successive upper-bound minimization (BSUM) method [29] which can naturally exploit the special structure of the unit modulus constraints. The PDD method has guaranteed convergence to KKT solutions of the spectral efficiency optimization problem. Moreover, it can be straightforwardly extended to the finite resolution phase shifter case. Our simulation results show that, even when both the transmitter and the receivers are equipped with the fewest RF chains that are required to support multi-stream transmission, hybrid precoding can still approach the performance of fully-digital precoding in both the infinite resolution phase shifter case and the finite resolution phase shifter case with several bits quantization.

Notations: scalars are denoted by lower-case letters, bold-face lower-case letters are used for vectors, and bold-face uppercase letters for matrices. For a matrix A, A T , A H , A † , Tr(A) and det(A) denote its transpose, conjugate transpose, pseudoinverse, trace, and determinant, respectively. I denotes an identity matrix whose dimension will be clear from the context. |x|, e {x} and x * are the absolute value, real part and conjugate of a complex scalar x, respectively, while x and X denote the Euclidean norm and the Frobenius norm of a complex vector x and a complex matrix X, respectively. x ∞ denotes the infinity norm. For a matrix X, we use [X] ij or X(i, j) to denote its (i, j)th entry. Particularly, we use X ∈ M to denote that each entry of X has a unit modulus constraint, i.e., |X ij | = 1. The distribution of a circularly symmetric complex Gaussian (CSCG) random vector variable with mean μ and covariance matrix C is denoted by CN (μ, C), and '∼' stands for 'distributed as'. C m ×n denotes the space of m × n complex matrices and R n denotes the n-dimensional real vector space.

Consider a narrow band single-cell mmWave downlink multiuser multi-stream MIMO system, where a base station (BS) equipped with N antennas and N RF (≤ N ) transmit RF chains sends signals to K > 1 users, each equipped with M > 1 antennas and M RF (≤ M ) receive RF chains. Let d > 1 denote the number of streams intended for each receiver. For successful symbol detection at each user, it is assumed that d ≤ min(N RF /K, M RF ). Furthermore, to reduce hardware complexity, we consider a hybrid digital and analog precoding architecture for both BS and users (see Fig. 1 of [15] for an illustration of the hybrid precoding architecture).

In the multi-user hybrid precoding scheme, the BS first processes the data streams digitally at the baseband using digital precoders, and then up-converts the digitally processed signals to the carrier frequency through RF chains, followed by an analog precoder which is implemented by analog phase shifters. Let V RF ∈ C N ×N R F denote the BS analog precoder, and V B B k ∈ C N R F ×d denote the digital precoder for user k's data stream s k ∈ C d×1 . Mathematically, the BS transmit signal after hybrid precoding is expressed as

and the received signal at user k is given by

where H k denotes the channel between the BS and user k, and n k denotes the additive white Gaussian noise (AWGN) with zero mean and variance σ 2 . At the receiver side, user k first processes the received signal by using an analog combiner U RF k ∈ C M ×M R F , and then down-converts the signals to the baseband through RF chains, followed by a digital combiner U B B k ∈ C M R F ×d to obtain the final processed signal, given by

(

Assuming Gaussian signaling for the data streams each with zero mean and unit variance, and that the noises n k 's and the data streams s k 's are independent of each other, the overall system spectrum efficiency maximization problem can be formulated as

) 1 The system rate expression does not include the digital combiners explicitly since it is well-known [30] that the optimal digital combiners (i.e., MMSE receivers) can achieve the maximum system rate. Moreover, for convenience, we use V to denote the set of variables V B B k 's and V R F ; similarly for U and other variables later. In addition, it is worth mentioning that, although rank constraints on precoders are meaningful, as in most previous works (e.g. [2], [6]- [9], [15], [30], [37], and [38]) we don't take them into consideration in our current formulation for simplicity and the complicated hybrid precoding problem with rank constraints on precoders is left open.

where the first constraint is the BS power constraint with power budget P ; the last two sets of unit modulus constraints are due to the fact that both the analog precoder and the analog combiners are implemented using low-cost phase shifters; and the matrix

is the so-called interference-plus-noise covariance matrix associated with user k.

It was shown in [15] and [18] that, when the number of transmit RF chains N RF are greater than or equal to twice the total number of streams, i.e., N RF ≥ 2Kd, the performance of the fully-digital precoding scheme can be perfectly achieved by the hybrid precoding scheme. This implies that problem (4) can be addressed by first solving the fully-digital precoding problem and then constructing the hybrid precoders V B B k and V RF based on the fully-digital precoder. It is noted that, the fully-digital precoding problem can be addressed using the well-known WMMSE algorithm [30], [31] and the corresponding hybrid precoders can be found in closed-form; see [15,Prop. 2]. However, for the general case, the hybrid precoding problem ( 4) is extremely hard to solve due to not only the coupling of the digital/analog precoders but also the unit modulus constraints. To the best of our knowledge, all the existing works that deal with these difficulties are heuristic algorithms and thus inevitably incurs performance loss. Although most of heuristic algorithms for single-user MIMO scenarios perform well, it is still not known how close their performance is to the optimal one. Hence, it is important to devise an algorithm that can solve problem (4) with some theoretical guarantees (e.g., achieve at least stationary solutions to problem (4)). In this paper, we propose using penalty dual decomposition method [28] to address problem (4). Before proceeding to solving problem (4), let us first give a brief introduction to the PDD method in what follows.

The PDD method is a general algorithmic framework that can be applied to the minimization of a nonconvex nonsmooth function subject to nonconvex coupling constraints. Considering that problem (4) has a differentiable objective function, we are here only concern about the differentiable case for ease of understanding of the PDD framework.

Consider the following problem

where f (x) is a scalar continuously differentiable function and h(x) ∈ R p×1 is a vector of p continuously differentiable functions; the feasible set X is the Cartesian product of n closed sets:

and n i=1 m i = m, and accordingly the optimization vari- (6) able x ∈ R m can be decomposed as

If no coupling constraints h(x) = 0 exist, the classical block coordination descent (BCD)-type algorithms [32] can be applied to decompose problem (P ) into a sequence of smallscale problems. This observation motivates us to dualize the difficult coupling constraints with appropriate penalty, and use coordinate-decomposition to perform fast computation, hence the name penalty dual decomposition method. Specifically, the PDD method applied to problem (P ) is a double-loop iterative algorithm. It employs inner iterations to solve to some accuracy a nonconvex augmented Lagrangian problem via an inexact or exact block coordinate descent method, while updating dual variables and a penalty parameter in outer iterations.

The PDD method is summarized in Table I, where the oracle 'BSUM(P k ,λ k , Lk , z k -1 , k )' means that, starting from z k -1 , the BSUM algorithm [29] or its variant randomized BSUM proposed in [28] for nonconvex constraint cases is invoked to iteratively solve problem (P k ,λ k ), given below

where L k (x) is the augmented Lagrange function with dual variable λ k and penalty parameter k . Further, the BSUM algorithm utilizes a locally tight upper bound of L k (x), denoted as Lk , and it terminates when certain solution accuracy k is reached. A basic implementation of the BSUM oracle is presented in Appendix A. Note that the detailed implementation and the convergence theory of the BSUM algorithm have been elaborated in [29]. Hence, although the BSUM oracle is the heart of the PDD method, we here omit its details for brevity but we will elaborate the BSUM algorithm when it is applied to the hybrid precoding subproblems later.

It is worth noting that, when the penalty parameter k in ( 7) is sufficiently small, the term h(x) 2 or h(x) ∞ will be driven to zero after solving problem (7). That is, in this case, solving problem (7) yields a solution to problem (P ) satisfy-ing the constraint h(x) ∞ = 0. However, we have generally h(x) ∞ = 0 during iterations. To measure the violation of the constraint h(x) = 0, we refer to the value of h(x) ∞ as constraint violation.

Define g(x) (g i (x i )) i . Then we have the following theory regarding the convergence the PDD method, which is a direct result of [28,Th. 3.1].

Corollary 3.1: Let {x k , ν k } be the sequence generated by Algorithm 1 for problem (P ), where ν k = (ν k i ) i denotes the Lagrange multipliers associated with the constraints g i (x i ) ≤ 0, ∀i. The stop criterion for the BSUM algorithm involved in Algorithm 1 is

with k , η k → 0 as k → ∞. Suppose that x is a limit point of the sequence {x k } and Robinson's condition holds for problem (P ) at x . Then x is a KKT point of problem (P ), i.e., it satisfies the KKT condition of problem (P ).

A general proof of this theorem can be found in [28]. Note that Robinson's condition is a constraint qualification condition which is often used to characterize the first order optimality condition of nonconvex problems [33]- [37]. For the problem at hand, it is equivalent to the commonly used Mangasarian-Fromovitz constraint qualification (MFCQ) condition. Further, a simple way to set η k is to make it explicitly related to the constraint violation of the last iteration or the current minimum constraint violation. For example, we set η k = 0.9 h(z k -1 ) ∞ in our simulations later. Furthermore, it is reasonable to terminate the BSUM algorithm based on the progress of the objective value

Here, the superscript 'r' denotes the BSUM iterations. In addition, since the penalty term h(x) ∞ vanishes eventually, a practical choice of the termination condition for the PDD method is h

Here, O is some prescribed small constant, e.g., O = 1e -6. See more details in [28].

This section proposes a PDD-based hybrid precoding method for spectral efficiency optimization. To tackle the difficulty arising from the rate function, i.e., the log det(•) function, we apply the PDD method to an equivalent problem of (4) instead of directly to (4).

Let us first rewrite problem (4) in the form of (P ). To make the power constraint easy to handle and remove the coupling of V RF and V B B k later through penalty [cf. eq. ( 13)], we introduce an auxiliary variable

Then we can recast problem (4) as max

Furthermore, to address the difficulty arising from the log det(•) function, based on the theory of the well-known WMMSE method [30], we write (10) in an equivalent WMMSE form, which is shown in Proposition 4.1 with the definition of meansquare-error (MSE) matrix

Moreover, if {U, V, W, X} is a KKT point of problem (12), {U, V, X} is a KKT point of problem (10), and further {U, V} is a KKT point of problem (4).

Proof: The proposition can be proven by following [30,Th. 3] and [38,Lemma 4.1]. We omit the proof for brevity.

Problem ( 12) is now in the form of problem (P ) is more tractable than problem (10) because the optimization with respect to W k 's is easy (see (17) below). Moreover, according to the result of Proposition 4.1, it is known that, the KKT solutions of problem (10) can be obtained by solving problem (12). Hence, in the following we solve problem (12) by using the PDD method [27].

Now we turn our attention to solving problem (12). It is readily known that the key to the PDD method is the inner iterations for solving augmented Lagrangian problems. Hence, our main efforts are devoted to developing BSUM/BCD algorithm for the augmented Lagrangian problem associated with (12). Specifically, by introducing dual variables Y k 's for the coupling constraints

It is seen that the constraints of the above problem are separable. Hence, we use BCD-type algorithm to address (13)

where

Since each matrix W k is positive definite, the optimal solution U B B k to problem (14) takes the form

The reason why we use pseudo-inverse here is that the matrix U H RF k A k U RF k may be rank-deficient during iterations.

2) The Subproblem w.r.t. W k : Consequently, given the optimal U B B k in (16), the optimal W k can be expressed as

where the last equality is obtained by plugging ( 16) into E k (U, X). Moreover, it can be verified that the bracketed matrix in the last equality is positive definite and so is the matrix W k .

3) The Subproblem w.r.t. V B B k : The variable V B B k appears only in the summand of the second term of the objective of (13), i.e., X k -V RF V B B k + ρY k 2 . Thus, the subproblem with respect to V B B k is given by

where Z k X k + ρY k . Solving the above unconstrained quadratic optimization problem, we obtain the optimal V B B k as follows

4) The Subproblem w.r.t. X: The variables X k 's appear in both the first and the second terms of the objective of (13). Thus, the subproblem with respect to X is given by

For notational convenience, we define

Then, with some appropriate rearrange, problem (20) can be equivalently written as

Further, by introducing a Lagrange multiplier μ ≥ 0 to the power constraint, we can express the optimal X k as

where the optimal multiplier μ can be easily found by using Bisection method such that

ii , U ρ is a unitary matrix consisting of the eigenvectors of A ρ and a i 's are the corresponding eigenvalues, i.e., U ρ diag {a 1 , a 2 , . . . , a N } U H ρ is the eigenvalue decomposition of A ρ .

5) The Subproblem w.r.t. V RF : The variable V RF only appears in the summand of the second term of the objective of (13). Thus, the subproblem with respect to V RF is given by (13) By appropriate rearrange, the above problem can be equivalently formulated as follows

where

Since the unit modulus constraints are separable, we can use one-iteration BCD-type algorithm to recursively solve problem (24). See the details in Appendix B.

6) The Subproblem w.r.t. U RF k : The variable U RF k appears only in the term E k (U, X). Thus, the subproblem with respect to U RF k is given by

By appropriate rearrange, the above problem can be equivalently formulated as follows

where A k is defined in (15),

Again, we can use one-iteration BCDtype algorithm to address problem (26). See the details in Appendix B.

The BCD-type algorithm for ( 13) is summarized in Table II. Note that, all updates can be implemented in parallel for different users. For instance, once U RF k 's and X k 's are given, each U B B k can be updated independently. By using parallel implementation, the proposed BCD-type algorithm is scalable to the number of users. In particular, when N >> M, it can be shown that the complexity of Algorithm 2 is dominated by Step 7, which is O(N 3 KI 2 ). Here, I 2 denotes the maximum number of iterations required by Algorithm 2.

In each iteration of the PDD method, after running Algorithm 2, we update the dual variable Y k or the penalty parameter ρ according to the constraint violation condition

, ∀k; otherwise, we increase the penalty parameter by updating ρ = cρ where 0 < c < 1. The PDD method is terminated when the feasibility of the penalized constraints (i.e.,

is approximately achieved; see the discussion in the end of Section III. In addition, it can be shown 2 that the MFCQ condition is satisfied for problem (12) at any feasible solution ({D X k }, D V R F , D U R F ). Hence, we can conclude that the PDD method converges to the set of KKT solutions of problem (4) according to the results of Proposition 4.1 and Corollary 3.1.

The PDD method requires running Algorithm 2 repeatedly. In this section, we provide a lightweight solution based on matrix approximation (MAP).

The intuition behind MAP is that [6], the performance achieved by hybrid precoding is upper bounded by the one achieved by the fully-digital precoding. Hence, if the optimal fully-digital precoder/decoder can be approximated well using hybrid precoder/decoder, we can obtain a good spectral efficiency performance that is close to that of the fully-digital precoding. Mathematically, given the optimal fully-digital precoder V opt k , ∀k, MAP for precoder design 3 is to find a structured matrix V k = V RF V B B k to approximate V opt,k , ∀k, i.e., solving the following matrix approximation problem

where

] and V RF ∈ M means that each entry of V RF has unit modulus. Note that, for tractability, we here have neglected the coupling power constraint which will be taken into consideration in the end by scaling V B B k 's to satisfy the power constraint. However, even after removing this difficulty, the problem is still difficult to solve due to the unit modulus constraint. The quadratic nature of the objective function motivates us to use BCD-type algorithm to address problem (27). Specifically, we divide the variable (V RF , V B B ) into N RF d + 1 blocks, i.e., V B B , V RF (i, j), i = 1, 2, . . . , N RF , j = 1, 2, . . . , d, and successively update each block by solving (27) for the selected block while fixing the others. Clearly, given V RF , V B B can be updated in closed-form:

Further, according to Appendix B, each entry of V RF , V RF (i, j), can be also updated in closed-form with the other variables fixed. Note that the proposed algorithm for problem (27) differs significantly from the alternating optimization method proposed in [7], where the analog precoder V RF is updated as a whole by using manifold optimization (MO) method [39]. Since the MO method is a gradient-descent-like iterative algorithm, our algorithm is more efficient than the alternating optimization method in [7] by using closed-form update of V RF . This will be verified later using numerical examples (see Fig. 3). 2 The MFCQ for this problem is equivalent to verify that 1) the equality constraint gradients are linearly independent, and 2) there exists

Note that both are easy to be verified. 3 The decoder can be designed in a similar way.

To summarize, the MAP method is a two-phase hybrid precoding method. In the first phase, we obtain the fully-digital precoder/decoder using the well-known WMMSE algorithm. In the second phase, we run BCD method to find the hybrid precoders/decoders. Assuming the number of transmit antennas N is much larger than the number of user receive antennas M , it can be shown that the complexity of the first phase is O(N 3 KI w m m se ) where I w m m se denotes the number of iterations required by the WMMSE method. Furthermore, as shown in Appendix B, updating all the entries of V RF once requires complexity of O(N 2 N 2 RF ). Similarly, updating U RF k 's requires complexity of O(KM 2 M 2 RF ). Hence, the complexity of the second phase is

where I bcd denotes the maximum number of iterations required for updating precoders/decoders. Therefore, the total complexity of the MAP method is cubic in the number of BS antennas. It is lower than that of the PDD method because the latter requires repeatedly running Algorithm 2 which also has cubic complexity in the number of BS antennas.

So far, we have assumed that arbitrary resolution phase shifters are available for realizing analog precoders. However, the hardware for accurate phase control could be very expensive. Furthermore, infinite resolution phase shifter is not always practical for large-array systems. Hence, we need to consider hybrid precoding with finite resolution phase shifters, i.e., the candidate phases for each phase shifter is finite.

Let F denote the set of finite phases, with |F| = 2 b , where b is the number of bits used to quantize the phases. The corresponding spectral efficiency optimization problem turns out to be max

Due to the combinatorial nature of the phases available for the phase shifters, problem (28) generally requires exhaustive search which however is computationally prohibitive in practice.

A heuristic way to deal with the constraints of finite phases is to first addressing the spectral efficiency optimization problem under the assumption of infinite resolution phase and then quantize each component of the analog precoders to the nearest point in the set F. This heuristic method could work effectively when the number of available phases is large (i.e., relatively high resolution phase). However, for the low resolution case (e.g., b ≤ 4), the effect of phase quantization could be significant. Hence, it is important to devise algorithms that can incorporate the constraints of finite resolution phases directly into the optimization procedure. Fortunately, our PDD algorithm 4 proposed above can be easily adapted to the finite resolution phase case. The modification lies only in Step 4 of Algorithm 4, i.e., modify Step 4 as follows x = arg min X(i,j )∈F e {b * X(i, j)}

The above problem can be globally solved via one-dimensional exhaustive search. As a result, the modified Algorithm 4 requires complexity of O(I 3 mn(mn + 2 b )).

Finally, we make a remark on iterations of our algorithms. All the algorithms proposed for both the infinite resolution phase case and the finite resolution phase case are iterative algorithms. While the algorithms requires a number of iterations for achieving convergence, we can terminate them early at the price of system spectral efficiency performance. It is worth mentioning that, when the algorithms are terminated early, we need to scale V B B k 's such the BS power constraint and eventually obtain a feasible solution to the spectral efficiency optimization problem.

This section presents simulation results to illustrate the performance of the proposed hybrid precoding algorithms. Since we focus on hybrid precoding for multi-stream multi-user MIMO cases without exploiting the angle information of channel knowledge, the simulation results for SU-MIMO, multi-user MISO, and single-stream multi-user MIMO cases are omitted in this paper and reported in [41]. The recent work [18] has shown that the regularized block diagonalization (RBD)-based hybrid precoding method performs better than ZF-based, BD-based, and MMSE-based hybrid precoding methods for multi-user MIMO cases. Hence, in the simulations, we compare our hybrid precoding methods with the RBD method in addition to the benchmark performance-the fully digital precoding schemes in terms of the achieved spectral efficiency. 5As in [15], we use a geometric channel model of L = 15 paths with uniform linear array antenna configurations and isotropic scattering. 6 Specifically, the channel matrix between the BS and each user is expressed as

where α k ∼ CN (0, 1) is the complex gain of the -th path, a r (φ k ) and a t (ϕ k ) are respectively the normalized receive and transmit array response vector at the azimuth angle of φ k ∈ [0, 2π) and ϕ k ∈ [0, 2π), which are given by

1, e j π sin(θ ) , . . . , e j (N -1)π sin(θ ) T . (30) In our simulations, it is assumed that the base station is equipped with N = 64 antennas and each user with M = 16 antennas. Unless otherwise specified, we set N RF = 8 and M RF = 4. For the PDD method, the initial penalty parameter ρ is set to 100/N and the control parameter c is set to 0.8. Furthermore, we set η 0 = 0 = 1e -3 and k = c k -1 . Moreover, in practical implementation, we set the maximum number of inner BSUM/BCD of the PDD methodfoot_3 to 30. The simulation results versus SNR are averaged over 100 channel realizations, where SNR is defined by SN R = 10 log 10 ( P σ 2 ).

In the first set of simulations, we examine the convergence performance of the proposed methods. To obtain a benchmark performance, we set K = d = 2 so that we have N RF = 2Kd and M RF = 2d. For this setup, it is known that the fully-digital precoding performance can be perfectly achieved by the hybrid precoding in the infinite resolution phase shifter case [15].

First, we examine the convergence performance of the PDD method. For 100 randomly generated multi-user channels, we run the PDD method for each channel with different random initialization. Note that, different problem instances yield different spectral efficiency value. To remove this effect and also clearly demonstrate how much percentage of fully-digital precoding performance the hybrid precoding can achieve, we normalize the objective value of problem (12) by the spectral efficiency value of the fully-digital precoding, and plot them in the left subfigure of Figs. 1 and2, where the average convergence behavior of the PDD method is shown for the case of infinite resolution phase shifters (denoted by b = ∞ for short) and finite resolution phase shifters with b = 1, respectively. In terms of the objective value and the constraint violation, it is observed from the plots that the PDD method can converge well within 50 iterations for both cases. Furthermore, it is seen that the hybrid precoding can achieve the same performance as the fully-digital precoding in the example of infinite resolution phase shifter case, implying its excellent convergence performance in achieving possibly optimal solutions. While for the finite resolution phase shifter case with b = 1, the hybrid precoding can achieve about 73.7 percent of the fully-digital precoding performance in this example.

Second, we examine the convergence performance of the MAP method as compared to the alternating optimization method in [7]. Since the two methods have the same way in the update of V B B but totally different ways in the update of V RF , we here focus on the comparison of the efficiency of the update of V RF and thus first compare the convergence performance of the BCD method and the MO method when they are applied to problem (27) with fixed V B B . Fig. 3 illustrates an average convergence behavior of the two methods over ten problem instances, where 'Approximation error' denotes the objective value of problem (27). It is seen that the BCD method requires significantly less iterations (generally serveral iterations) for convergence than the MO method, while both can achieve the same convergence result. Furthermore, Fig. 4 shows that, with different initializations, the MAP method can always achieve global optimality of problem (27). This will be further verified by the simulation results later.

In the second set of simulations, we demonstrate the spectral efficiency performance of the proposed hybrid precoding methods as compared to the benchmark fully-digital precoding performance and the existing multi-user hybrid precoding method-RBD [18] in the infinite resolution phase shifter case. Fig. 5 shows the spectral efficiency performance of various precoding methods for the case of K = d = 2. In this case, it is known that the fully-digital precoding performance is achievable by the hybrid precoding. Fig. 5 shows that both the PDD method and the MAP method can achieve the same performance as the fully-digital precoding. This again verifies the excellent convergence performance of the PDD method and the MAP method (possibly achieve global optimality in this case). Moreover, it is observed from Fig. 5 that both the PDD method and the MAP method has better performance than the RBD method with about 1 dB gain.

Fig. 6 shows the spectral efficiency performance of various hybrid precoding methods for the cases when N RF = Kd < 2Kd. It is seen that, even for the case when N RF < 2Kd, the PDD method can achieve a performance that is extremely close to the FD precoding performance. Furthermore, one can see that, unlike the case of N RF ≥ 2Kd (where the MAP method has the same performance as the PDD method), the PDD method has a better performance than the MAP method in the cases of N RF < 2Kd. Moreover, it is observed again that both the PDD method and the MAP method outperform the RBD method, and the performance gaps increase with the SNR. The reason for this observation is explained as follows. Recall that both the MAP method and the RBD method first neglects the power constraint and then scale the digital precoder V B B k 's to satisfy the power constraint. Such a heuristic scaling approach inevitably incurs performance degradation especially when P is large (i.e., the SNR is large). Particularly, when SN R = 6, the PDD method improves the spectral efficiency of the MAP method and the RBD method by 5 bps/Hz and 8 bps/Hz, respectively.

In the third set of simulations, we demonstrate the spectral efficiency performance achieved by the PDD-based hybrid precoding method in the finite resolution phase shifter case. Since it is not possible to get good matrix approximations of the fullydigital precoder/decoders in this case, the MAP method does not work well and thus its simulation result is not shown here. Furthermore, since the RBD method was proposed only for the infinite resolution phase shifter case in [18], we adapt it to the finite resolution phase shifter case by quantizing the analog precoder/decodersfoot_4 followed by the design of the digital precoder/decorders. For convenience, we still refer to the modified RBD method in the finite resolution phase shifter case as 'RBD' in the plot. In addition, since the hybrid precoding problem has much more local maxima in the finite resolution phase shifter case than in the infinite resolution phase shifter case, the PDD method could get easily trapped in some bad point and thus may result in bad performance. To deal with this issue, we run 20 iterations of the PDD method assuming infinite resolution phase shifters to get a good initialization, followed by the PDD method of the finite resolution phase shifter case. 9 Such a strategy makes the PDD method work very well in the finite resolution phase shifter case from our numerical experience. Fig. 7 illustrates the spectral efficiency of the PDD method and the modified RBD method when three bits quantization is used, i.e., b = 3. It is observed that the PDD method can still achieve most of the spectral efficiency of the FD precoding in this example of finite resolution phase shifter case. Moreover, one can see again that the PDD method outperforms the RBD method and their gap increases with the SNR. For example, the gap between the PDD method and the modified RBD method is about 5 bps/Hz when SN R = 6 dB, while it is much smaller (almost negligible) when SN R = -10 dB.

Fig. 8 illustrates the spectral efficiency of the PDD method versus SNR when different quantization levels are used. It is observed that, the spectral efficiency of the PDD method improves when more bits are used in the phase quantization, and the improvement shrinks as the quantization level increases. Particularly, one can see that three bits quantization is enough for achieving about 95% of the fully-digital precoding performance in this example.

Fig. 9 shows the spectral efficiency of the PDD method when different number of RF chains are used with the lowest resolution phase shifters, i.e., b = 1. It is seen that the performance gap between the PDD-based hybrid precoding and the FD precoding can be reduced by increasing the number of RF chains. Therefore, the number of RF chains can be used to trade off the resolution of the phase shifters in hybrid precoding design.

Lastly, we consider the cases when the numbers of transmit/receive RF chains fullfil the minimum requirement, i.e., N RF = Kd and M RF = d, and examine how much percentage of the FD precoding performance the hybrid precoding can achieve. In the simulations, we test five system setups with different combination of K and d which satisfies N RF = Kd and M RF = d. For each setup, we run the PDD method for 100 randomly generated channels under three quantization levels b = ∞, b = 4, and b = 2. The minimum, average, maximum spectral efficiency of the PDD method relative to the spectral efficiency of the FD precoding are respectively listed in Table I for SN R = 0 dB. It can be observed that, in the infinite resolution phase shifter case, the hybrid precoding can achieve more than 95 percentage of the FD precoding performance, implying the excellent efficiency of the hybrid precoding while enjoying the benefit of significantly reducing the number of RF chains in Massive MIMO systems. Furthermore, even with low resolution phase shifters, the hybrid precoding can still achieve most percentage of the FD precoding performance, e.g., it is on average about 95% when b = 4 and about 80% when b = 2.

By applying penalty dual decomposition method, this paper has proposed an iterative algorithm to address the difficult hybrid precoding problem for mmWave multi-user MIMO systems. Different from the existing hybrid precoding algorithms which are all heuristic, the proposed algorithm has guaranteed convergence to KKT solutions of the hybrid precoding problem. Although the optimality of the proposed hybrid precoding method is not proven, simulation results verify that the proposed hybrid precoding method performs very well and its achievable spectral efficiency is very close to the performance of the fully-digital precoding, implying the capability of hybrid precoding in achieving near-optimal spectral efficiency while greatly reducing the number of RF chains. It is worth mentioning that, all the results are obtained under the assumption of perfect channel state information. In the future, we will consider hybrid precoding with imperfect channel state information.

We outline the basic BSUM algorithm [29] in Table IV which can be used to implement the oracle x k = BSUM(P k ,λ k , Lk , x k -1 , k ) in Step 2 of the PDD method. Here, Lk (y i ; w) denotes a locally tight upper bound of L k (x) w.r.t x i (i.e., the i-th block of x) at the point w.

In the BSUM algorithm [29] applied to a minimization problem with multiple block variables, each time one block variable is cyclicly picked to be optimized while fixing the others by minimizing a locally tight upper bound of the objective. In particular, when the upper bound is simply chosen as the objective function itself, the BSUM algorithm reduces to the BCD method [32]. Hence, the former includes the latter as a special case. On the other hand, note that the BSUM algorithm [29] was proposed for convex constraint cases. When the problem has separable nonconvex constraints, we need to resort to a randomized BSUM algorithm proposed in [28] instead of BSUM for theoretical convergence guarantee. For convenience, we sill refer to it as BSUM algorithm, but keep in mind that the BCD or BSUM algorithm used throughout this paper refers to randomized BSUM algorithm in nonconvex constraint cases, where each time one block variable is randomly picked to be optimized while fixing the others.

This appendix provides an iterative algorithm to address the following quadratic optimization problem with unit modulus (31) constraints: min X∈M φ(X) Tr(X H AXC) -2 e Tr(X H B) (31) where the matrices A ∈ C m ×m and C ∈ C n ×n are positive semidefinite and M denotes the set of unit modulus constraints, i.e., |X(i, j)| = 1, ∀i, j.

We use BCD-type algorithm to address problem (31) with guaranteed convergence to stationary solutions [32], i.e., in each step we update one entry of X while fixing the others. Without loss of generality, let us consider the problem of minimizing φ(X) with respect to X(i, j) subject to the unit modulus constraint |X(i, j)| = 1, i.e., min |X(i,j )|=1 φ(X) (32) It is easily known that, the function φ(X) with respect to X(i, j) can be expressed as a quadratic function of X(i, j) in the form of φ(X(i, j)) a|X(i, j)| 2 -2 e {b * X(i, j)} for some real number a and some complex number b. Considering that |X(i, j)| = 1, problem (32) reduces to max |X(i,j )|=1 e {b * X(i, j)} . (33) It follows that the optimal X(i, j) is equal to b/|b|. Hence, in order to update X(i, j), we only need to know the value of b.

In what follows, we show how the complex number b can be easily obtained. First, we have [40] ∂ φ(X(i, j)) ∂X * (i, j) X(i,j )= X(i,j ) = 1 2 (a X(i, j)b).

On the other hand, we have [40] ∂φ(X) ∂X * X= X = 1 2 (A -B).

Combining the above equations, we obtain [A XC -B] ij = a X(i, j)b. By expanding [A XC] ij and checking the coefficient of X(i, j), we have a X(i, j) = A(i, i) X(i, j)C(j, j). It follows that b = A(i, i) X(i, j)C(j, j) -[A XC] ij + B(i, j).

According to the above analysis, the entries of X can be recursively updated. The corresponding algorithm for updating X is summarized in Table V, where the recursion step 5 is due to the fact that Q should be updated accordingly once X(i, j) is updated (which is done in step 6). It is easily known that step 5 is the most costly step requiring complexity O(mn). Hence, it can be shown that Algorithm 4 has complexity of O(I 3 m 2 n 2 ) where I 3 denotes the total number of iterations required by Algorithm 4.