1

Joint User Scheduling and Beam Selection
Optimization for Beam-Based Massive MIMO
Downlinks
arXiv:1711.06915v1 [] 18 Nov 2017

Zhiyuan Jiang, Member, IEEE, Sheng Chen, Sheng Zhou, Member, IEEE,
Zhisheng Niu, Fellow, IEEE

Abstract
In beam-based massive multiple-input multiple-output systems, signals are processed spatially in the
radio-frequency (RF) front-end and thereby the number of RF chains can be reduced to save hardware
cost, power consumptions and pilot overhead. Most existing work focuses on how to select, or design
analog beams to achieve performance close to full digital systems. However, since beams are strongly
correlated (directed) to certain users, the selection of beams and scheduling of users should be jointly
considered. In this paper, we formulate the joint user scheduling and beam selection problem based on
the Lyapunov-drift optimization framework and obtain the optimal scheduling policy in a closed-form.
For reduced overhead and computational cost, the proposed scheduling schemes are based only upon
statistical channel state information. Towards this end asymptotic expressions of the downlink broadcast
channel capacity are derived. To address the weighted sum rate maximization problem in the Lyapunov
optimization, an algorithm based on block coordinated update is proposed and proved to converge to the
optimum of the relaxed problem. To further reduce the complexity, an incremental greedy scheduling
algorithm is also proposed, whose performance is proved to be bounded within a constant multiplicative
factor. Simulation results based on widely-used spatial channel models are given. It is shown that the
proposed schemes are close to optimal, and outperform several state-of-the-art schemes.

Index Terms
Massive MIMO, user-scheduling, hybrid beamforming, statistical CSI
Z. Jiang, S. Chen, S. Zhou and Z. Niu are with Tsinghua National Laboratory for Information Science and Technology,
Tsinghua University, Beijing 100084, China. Emails: {zhiyuan@, chen-s16@mails., sheng.zhou@, niuzhs@}tsinghua.edu.cn.
This work is sponsored in part by the Nature Science Foundation of China (No. 61701275, No. 91638204, No. 61571265,
No. 61621091), the China Postdoctoral Science Foundation, and Intel Collaborative Research Institute for Mobile Networking
and Computing. The corresponding author is Dr. Sheng Zhou.

2

I. I NTRODUCTION
In massive multiple-input multiple-output (MIMO) based wireless communication systems,
the spectral and radiated energy efficiency can be both boosted by the deployment of massive
number of antennas [1]. Moreover, the high beamforming gain of a massive antenna array is
the main enabler for millimeter-wave systems against high pathloss. Therefore, it is extremely
important to design high-performance, efficient and practical transmission strategy in massive
MIMO systems for the emerging 5G cellular system.
Under the assumption that full digital signal processing is performed at the base station (BS)
side with massive antenna arrays, the system performance has been widely investigated, e.g.,
in [1]–[3]. However, it is widely accepted that full digital signal processing implementation
encounters very severe challenges in practice, on account of the following impediments.
Radio-frequency (RF) chain hardware cost and power consumptions. Full digital signal
processing requires that all antennas can be digitally controlled from baseband. Hence, one
dedicated RF chain, including e.g., low-noise amplifier, analog-digital-converter (ADC), power
amplifier and etc., is needed for each antenna. In massive MIMO systems, not only is this
requirement entails a dramatic increase in the deployment cost of the system, but also that the
power consumption would be driven up to a prohibitive level. As indicated in the previous work
[4] [5], concretely, a BS with 256 RF chains consumes about 10 times the power (only the RF
chains) as compared with an entire current long-term-evolution (LTE) BS.
Baseband signal processing complexity. The spatial baseband processing includes multiple
kinds of matrix operations, such as inversions and singular-value-decompositions (SVDs) whose
complexity scales with M 3 where M is the number of antenna elements for full digital processing. Moreover, these extremely demanding matrix operations are required to be executed very
frequently (once every 1 ms for spatial precoding in LTE systems). This is very challenging to
the design of baseband processing units, both in terms of chip costs and power consumption.
System specific limitations. Aside from the first two challenges, there are some other practical
considerations which are system-specific. For example, the fronthaul interface in cloud radio
access networks (C-RAN) poses a serious limitation in the number of data streams that can be
transmitted between the remote-radio-units (RRUs) and the baseband units (BBUs). Considerable
amounts of work has been dedicated to the signal spatial compression in C-RAN [6]. Moreover,
the channel state information (CSI) acquisition overhead in frequency-division-duplexing (FDD)

3

RF Chain 1

UE 1

Ant 1

Queue 1
RF Chain 2

Queue 2

User-Beam
Mapping
(User
scheduling)

…

UE 2

Antenna-Beam
Mapping
(Beam
selection)

Ant 2

RF Chain K

UE 3

Queue 3
User and beam
correspondence

Queue N

…

…
UE N

Ant 3

Ant M

Fig. 1. An illustration of the beam-based massive MIMO systems where user scheduling and beam selection are correlated.

system scales with the number of digitally controllable antennas. It constitutes a major bottleneck
in realizing the massive MIMO gain in FDD systems.
In view of these challenges, architectures with low RF- and processing-complexity have
been proposed extensively, e.g., in [7]–[13]. The existing literature can be divided into three
categories. The first is hybrid beamforming, which adopts an RF front end with an analog
beamforming module such that the number of RF chains is significantly reduced [8]. Although
the analog beamforming module is usually composed of phase shifters with constant-amplitude
beamforming weights to save hardware cost, the high-speed phase shifters, whose quantity is the
same with the number of antenna elements, pose a drastic challenge to the cost of RF front ends.
In this regard, the recently proposed beamspace MIMO architecture [10] adopts a lens antenna
array which acts analogously like a lens focusing on light beams from different directions. It
transforms the signal into the angular domain and thus reduces the number of RF chains due
to signal angular sparsity. Since it does not require any phase shifters, the total cost is reduced,
and therefore it is considered to be one of the candidate solutions to the 5G millimeter-wave
massive MIMO systems. The other approach is based on digital beamforming which involves
multi-layer signal processing [4], [11], [14]. Although the number of RF chains is not reduced,
the processing complexity and pilot overhead problems are partly solved.
In essence, all the aforementioned solutions aim at providing comparable performance as
full digital processing systems with limited number of RF chains and reduced complexity in
massive MIMO systems. Based on Fig. 1, the current literature mainly focuses on the right side

4

of the figure, i.e., the antenna to beam mapping and beam selection schemes, which leverages
the angular domain power sparsity of the channel to transform the signals from the antenna
domain to the beam domain. Towards this end, beam sweeping and steering methods in the
hybrid beamforming architecture can be used to capture the signal direction [15]. Other methods
adopt joint analog and digital precoding design [9] [16]. In beamspace MIMO systems, the
lens antenna array can be regarded as a directional beamforming module with low cost. On the
other hand, the left side of the Fig. 1 which represents user-beam mapping, is scantly treated.
The user-beam mapping essentially deals with user scheduling in the beam domain. Unlike
the previous user scheduling related work, e.g., in [17]–[19], the user scheduling problem in
the beam domain is tangled with the beam selection. In reality, due to the angular sparsity of
the massive MIMO channel [20], the beams, which represent the signal directions, are strongly
related to the users, in the sense that each beam usually contains signals of very few (possibly
one) users. Therefore, the user scheduling and beam selection have to be jointly considered to
avoid possible performance degradation due to user-beam mismatch.
The channel state information (CSI) is of vital importance to the system. The CSI can be
categorized as instantaneous CSI and statistical CSI. It is worthwhile to emphasize the specific
CSI usage at each stage (time scale) of the beam-based massive MIMO transmissions since most
existing work ignores this and assumes instantaneous CSI is always available [21] [22] . We
propose that beam-based downlink scheduling should be performed only based on the statistical
CSI. The reason is two-fold. From an implementation perspective, the statistical CSI is much
easier to obtain than instantaneous CSI, attributing to the fact that statistical CSI can be obtained
with much lower cost because a) it can be estimated without dedicated pilots [23]; b) it varies
at a lower speed (in the order of 1 second to 10 seconds) compared with instantaneous CSI
(in the order of 1 ms to 100 ms) [24]. Moreover, in C-RAN systems, the beamforming module
is integrated with the remote radio heads (RRHs) and hence limited computation capability is
expected [25] which prevents us from using complicated channel estimation schemes. On the
other hand, it is also theoretically possible to only rely on statistical CSI in the user scheduling
and beam selection phase, since beams are essentially long-term statistics. Furthermore, the
statistical CSI can be obtained efficiently with a limited number of RF chains based on, e.g.,
compressive sensing based channel estimation schemes [7] [26].
In this paper, we aim to address the user and beam joint scheduling problem in beam-based
massive MIMO downlinks. The contributions include:

5

•

We formulate the problem based on the Lyapunov-drift optimization framework. An optimal
scheduling policy is proposed thereby to achieve optimum utilities. The optimality proof is
given which shows the achieved utility is arbitrarily close to the optimum.

•

To address the queue weighted sum rate maximization problem arisen in optimizing the
Lyapunov-drift, which is a mixed integer programming (MIP) problem, the block-coordinatedupdate-based (BCU-based) algorithm which deals with the continuous convex relaxation of
the MIP problem is proposed. In order to implement the algorithm based on statistical CSI,
a deterministic equivalence of the downlink broadcast channel capacity in the large antenna
array regime is derived, depending only on statistical CSI. The BCU-based algorithm is
proved to converge to the global optimum of the relaxed problem. An iterative water-filling
based approach is also proposed to reduce the number of iterations.

•

Furthermore, a low-complexity incremental greedy algorithm is proposed. We prove that
the greedy algorithm can achieve near-optimal performance, within a multiplicative factor
due to the submodular property of the problem.

•

By simulations, it is shown that the proposed algorithms outperforms several state-of-theart beam selection schemes. Moreover, since it is based on statistical CSI, the frequency of
executing the algorithm is significantly reduced, making it more preferable in practice.

A. Related Work
The proposed joint user and scheduling schemes are related to the beam selection problem in
beamspace MIMO systems [21], [22], [27]–[30], or more generally antenna selection problem
in MIMO systems [31]. However, the considered joint scheduling problem in beam domain is
unique, in the sense that the beam magnitudes are strongly correlated with users. The beamuser scheduling problem is also considered in a switched-beam based massive MIMO system
in [28], where one pre-defined (fixed) beam is associated with one user. In [29], a greedy
joint scheduling of beams of users is proposed and we will compare our results with it in the
simulations. Concerning the literature related to the mathematical treatment adopted in the paper,
the Lyapunov-drift optimization framework is attributed to the pioneer work by Neely [32]. The
large system deterministic equivalence to derive the downlink channel capacity is related to the
celebrated random matrix theory [33]. Furthermore, in the algorithm design, the BCU technique
dates back to multi-convex optimization, e.g., in [34], and the approximation factor of the greedy
algorithm is related to the submodular set function optimization problem as in [35].

6

B. Paper Organizations and Notations
The remainder of the paper is organized as follow. In Section II, the system model and channel
model are presented, and the problem is formulated. In Section III, the Lyapunov-drift approach
is used to design an optimal scheduling policy. In Section IV, a BCU-based scheduling algorithm
is described to address the queue weighted sum rate maximization problem. In Section V, a lowcomplexity greedy algorithm is presented. The simulation results are conveyed in Section VI.
Finally, we conclude our work in Section VII.
Throughout the paper, we use boldface uppercase letters, boldface lowercase letters and
lowercase letters to designate matrices, column vectors and scalars, respectively. X T and X †
denotes the transpose and complex conjugate transpose of matrix X, respectively. Xi,j and xi
denotes the (i, j)-th entry and i-th element of matrix X and vector x, respectively. tr(X) denotes
the trace of matrix X. Denote by E(·) as the expectation operation. Denote by IN as the N
dimensional identity matrix. The logarithm log(x) denotes the binary logarithm.
II. S YSTEM M ODEL AND P ROBLEM F ORMULATION
A. Signal Model
The single-cell system downlink is considered in this paper, where a BS with M co-located
antennas transmits to Nt single-antenna users.1 The BS has K RF chains (K ≤ M ). Assuming
narrow-band and time-invariant channels2 , the receive signal of user-n is, ,
yn = h†n x + nn ,

(1)

where hn is an M -dimensional channel vector, x is the downlink transmit signals, and nn denotes
the i.i.d. Gaussian additive noise with unit variances. The downlink channel matrix is denoted
by H = [h1 , h2 , ..., hNt ]† . The transmit signal after beamforming can be written as
x = Ba s,
1

(2)

The proposed schemes can be straightforwardly extended to multiple-antenna-user case by treating multiple antennas of a

user as multiple users with an identical channel correlation matrix [23] and setting the backlog pressure Qn (t) in P3 to be the
same for these antennas.
2

Wideband channels can be decomposed to a set of parallel narrow-band channels by, e.g., orthogonal-frequency-division-

multiplexing (OFDM) modulations. The time-invariant channel assumption essentially deals with data transmission inside the
channel coherence time (or block length in the block fading channel model).

7

where s denotes the K-dimensional digitally precoded data symbols for the scheduled, i.e.,
spatial-multiplexed Ns users (Ns ≤ K such that a linear digital precoder such as zero-forcing
precoder can eliminate the inter-user-interference,3 and obviously Ns ≤ Nt ). The RF (analog)
beamforming, which can be realized by the lens antenna array and beam selection in beamspace
MIMO or general analog beamforming in hybrid beamforming architectures, is denoted by Ba
with dimension M × K. On account of the analog beamforming, the effective channel observed
from baseband is
H̄ = HBa ,

(3)


†
where the effective channel vector corresponding to user-n is denoted by h̄n and H̄ = h̄1 , ..., h̄Nt .
The RF beamforming considered in the paper is a widely adopted directional beamforming
scheme4 , and hence
Ba = BDFT Σb ,

(4)

where BDFT is the equivalent discrete-Fourier-transform (DFT) matrix (or Kronecker product
of DFT matrices for uniform planar antenna arrays (UPAs)). The beam selection decision is
denoted by the diagonal matrix Σb whose entries are binary, i.e., (Σb )i,i ∈ {0, 1}, ∀i.
B. Channel Model
Based on a geometry-based channel model [36], the channel vector of the n-th user can be
written as
r
hn =

Ln


M X
(n)
(n)
(n)
βl α θl , ψl
,
Ln l=1

(5)

where Ln denotes the total number of multi-path components (MPCs) in the propagation medium
including line-of-sight (LoS) and none-line-of-sight (NLoS) MPCs. The amplitude of each MPC
(n)

(n)

is denoted by βl , and θl

(n)

and ψl

denote the azimuth and elevation angles of the l-th arriving

3

In the simulations, we assume K = Ns and hence K ≤ Nt .

4

Note that the RF beamforming in this paper can be readily generalized to arbitrary beam pattern, e.g., the eigenvector-based

beam pattern in [14], by replacing the DFT-based beamforming matrix BDFT with the desired beamforming matrix. Also note that
although directional beamforming is adopted, the proposed schemes can adapt to the channel variation more flexibly compared
with traditional directional antenna based systems, by updating the channel statistics estimations and adjusting beam patterns.

8

MPC, respectively. Thereby, the steering vector for one MPC is given by (assuming UPA whereas
one-dimensional uniform-linear-array (ULA) can be regarded as a special case)





(n)
(n)
(n)
(n)
1
−j2π mλh sin θl cos ψl +nλv cos θl sin ψl
(n)
(n)
UPA
1, ..., e
, ...,
α
θl , ψl
=√
M

 T
(n)
(n)
(n)
(n)
−j2π (H−1)λh sin θl cos ψl +(V −1)λv cos θl sin ψl
e
,

(6)

where H and V denote the number of columns and rows in the UPA, respectively, and λh and
λv are the antenna spacing in the horizontal and vertical domains, respectively. The order of
elements in the steering vector is mapped to the indexing order of antennas in the UPA. The
physical meaning of the steering vector and channel representation in (5) is that for an MPC


(n)
(n)
with direction-of-arrival (DoA) θl , ψl , the array response is given by (6). Summing up
all the contributing MPCs, we obtain the compound channel representation in (5). Based on the
channel model, the RF beamforming in (4) can take advantage of the limited number of MPCs
as compared with the number of antennas, and only selects a subset of the beams to attain
equivalent performance (signal power) with a smaller number of RF chains.
However, the beam selection cannot be isolated from the user scheduling problem. Apart
from the reasons given in Section I, from a throughput perspective, different users have different
transmission needs resulting from traffic demand or fairness considerations. Therefore, the user
scheduling and beam-domain CSI should also be jointly considered. The joint beam-domain
massive MIMO downlink scheduling problem is formulated as follows.
C. Problem Formulations
The long-time average rate of user n is denoted by R̄n , and the instantaneous rate of user n
at time t is denoted by Rn (H(t), π(t)), given the channel coefficients H(t) and control policy
(user scheduling and beam selection as far as the paper is concerned) π(t). Note that this does
not mean the scheduling decision relies on the availability of instantaneous CSI, as in Proposition
1 an deterministic equivalence of the rate expression will be derived which is based solely upon
statistical CSI. Based on ergodicity, R̄n = E{Rn (H, π)}, ∀n ∈ {1, ..., Nt }, where the expectation
is taken over channel coefficients H(t) and possibly π(t) when a stochastic control policy is
considered. The achievable ergodic rate region can be characterized as
[
R̄ : 0 ≤ R̄n ≤ E [Rn (H, π)] ,
R = coh
π∈X

(7)

9

where R is a Nt -dimensional region, R̄n is its n-th component, and “coh” denotes the closure of a
convex hull. The set of all feasible scheduling policies is denoted by X . The utility maximization
problem is formulated as

maximize U R̄ , subject to R̄ ∈ R,

P1:

(8)

Σu ,Σb

∆

where Σu is a diagonal matrix denoting user scheduling decision at time t, i.e., si = (Σu )i,i ∈

∆
{0, 1}, and bi = (Σb )i,i denotes the beam selection decision. The network utility function U R̄
is defined as a function of the long-time average rate for each user, e.g.,


Usum R̄ =

Nt
X

R̄n

(9)

log(R̄n + cn )

(10)

n=1

for sum rate maximization,


Upfs R̄ =

Nt
X
n=1

for proportional-fairness scheduling (PFS) [37], where cn ’s are non-negative constants to regularize the logarithm expressions, and typical value is cn = 0, ∀n ∈ {1, ..., Nt } for exact PFS
or cn = 1, ∀n ∈ {1, ..., Nt } to ensure positive objective function value which is a mathematical

convenience. Basic properties of the utility function U R̄ are required, e.g., concavity and
monotonicity [38], over the rate vector (R1 , ..., RNt ).
III. O PTIMAL B EAM -BASED J OINT S CHEDULING P OLICY
To solve P1, it is found that two severe challenges exist. First, the ergodic capacity region
R does not yield a closed-form expression. The work in [39] [40] characterizes the broadcast
channel (BC) capacity region and a duality between BC and multiple-access-channel (MAC)
in the sense of both capacity region and outage probability is found. Moreover, an iterative
water-filling algorithm is proposed to calculate it, given the instantaneous channel coefficients.
Nevertheless, no closed-form expressions are available except for capacity bounds [41].5 Secondly, the scheduling and beam selection decisions should be made dynamically to match the
channel variations and user traffic in time. To address these issues, we seek to leverage a powerful
tool of Lyapunov-drift optimization which is shown to have superior performance compared to
5

The broadcast channel capacity formula is adopted as the optimization objective in this paper due to its better generality

compared with, e.g., achievable rates based on linear precoding schemes. It is also because that non-linear downlink transmission
schemes, e.g., non-orthogonal multiple access schemes, are attracting more and more attention recently.

10

static solutions [42] with simple decision structures (max-weight structure [43]); thereby, P1 is
decomposed into P2 and P3 described in the following.
Essentially, P1 is a time-average network utility maximization problem which is hard to
solve directly. The Lyapunov-drift approach decomposes the time-average optimization into
optimization in each scheduling step; and the resultant sub-problems are formulated as P2 and
P3. By applying the solutions to P2 and P3 at each scheduling step, the time-average network
utility can be optimized.
A. Lyapunov-Drift Based Network Utility Maximization
To maximize the network utility function in (8), the transmission need of each user, which
is determined by the transmission history and utility function, is represented by a set of virtual
queues. The arrival process is designed to reflect the transmission urgency of each user and
a max-weight algorithm is applied to stabilize the queues whenever possible. Specifically, let
Qn (t) denote the virtual queue length in bits of user n at the beginning of t-th scheduling step,
let an (t) denote the number of (virtual) arrival bits which are optimization variables for utility
maximization, and let µn (t) denote the allocated number of service bits to queue-n, which equals
the allocated number of service bits between scheduling steps. The queuing dynamics are written
as
Qn (t + 1) = Qn (t) − µ̃n (t) + an (t),
where µn (t) =

PT

τ =1

(11)

Rn (H(τ ), π(τ )), the number of channel uses between scheduling steps

is denoted by T , and µ̃n (t) = min{Qn (t), µn (t)} denotes the number of actual service bits,
considering the circumstances that sometimes the queue is emptied given the amount of allocated
service bits. Notice that the queues here are created virtually to facilitate the utility maximization
and thus they are not real traffic patterns. In Section VI (Fig. 5), we extend to stochastic real
traffic scenarios in simulations. The optimal beam-based downlink scheduling policy at a given
scheduling time t, i.e., a dynamic policy which achieves the solution to (8), can be described as
below.
Admission control: For virtual queue Q(t) = [Q1 (t), ..., QNt (t)], let the number of arrival
bits, i.e., a(t), be the solution of
P2:

maximize V U (a(t)) − a(t)T Q(t),
a(t)

subject to 0 ≤ an (t) ≤ Amax , ∀n ∈ {1, ..., Nt },
(12)

11

where V and Amax are pre-defined constants6 .
Scheduling: Given the arrival process determined above, the service, i.e., the joint scheduling
decisions, is based on the solution of the following problem:
#
"
Nt
T
X
X
Qn (t)sn
Rn (H(τ ), Σu , Σb , p)
P3:
maximize
Σu ,Σb ,p

s.t.,

n=1
Nt
X

(13)

τ =1

pn ≤ P,

(14)

n=1
Nt
X

si = Ns ,

i=1

M
X

bi = K, si ∈ {0, 1}, bi ∈ {0, 1},

(15)

i=1

where p = [p1 , ..., pNt ] denotes the transmit power corresponding to Nt user data streams and
hence P in (14) is the sum power constraint. The scheduling decisions are denoted by binary
variables si and bi . Σb and Σu are diagonal matrices consisting of si and bi , respectively. The
downlink instantaneous transmission rate Rn (H(t), Σu , Σb , p) is a function of the downlink
transmit power allocation, channel coefficients, and scheduling decisions. The departure from
P
the n-th virtual queue is µn (t) = sn Tτ=1 Rn (H(τ ), Σu , Σb , p).
It is observed that the admission control problem P2 is a convex problem and hence is easy
to solve. For instance with PFS utility, the optimal admission control is given by


V
∗
an (t) = min
, Amax , n ∈ {1, ..., Nt }.
Qn (t)

(16)

However, the problem P3 is an MIP problem, which is NP-complete [44]. Before diving into
details on solving P3 in the following sections, we assume the optimal solutions to both problems
are obtained for the moment, which is denoted by π ∗ . The optimality of the algorithm is
established in the following theorem.
Theorem 1: Denote

R̄∗ = arg max U R̄ .

(17)

R̄∈R

Suppose the transmission rate is bounded, i.e., 0 ≤ R̄n ≤ Rmax , ∀n ∈ {1, ..., Nt }, the utility
function U (·) is concave and entry-wise non-decreasing, and bounded on [0, Rmax ]. The channel
6

For Typical values, V and Amax can be approximately 100-fold of the service rate.

12

coefficients H(t) are i.i.d. over different scheduling periods, then based on the scheduling
algorithm resulting from P2 and P3, the following conditions are met.
!
τ −1

1X
E[R(t)] ≥ U R̄∗ − C/V,
lim inf U
τ →∞
τ t=0

(18)

E [Qn (τ )]
= 0, ∀n.
τ →∞
τ
lim

(19)

Proof: See Appendix A.
Remark 1: Theorem 1 reveals that the utility function of the time-averaged transmission rate
based on the scheduling decisions derived in P2 and P3 is within a constant (arbitrary small if
V is large) to the optimum and the virtual queues are mean-rate-stable, where C is a constant
related to Amax (40). In the following, we will dig into the methods to solve P3.
IV. B LOCK C OORDINATE U PDATE BASED M ETHOD FOR P3
This section is dedicated to solving the scheduling problem of P3 only based on the knowledge
of statistical CSI. The previous section establishes the optimality of the proposed beam-based
scheduling algorithm given the solutions of P2 (generally easy to solve) and P3. However, due
to the NP-hardness of P3, explicit solutions are hard to attain. More importantly, it is proposed
that the scheduling decisions of P3 should only rely on statistical CSI, rendering the solution
even more intractable. Towards this end, an algorithm based on solving the convex relaxation
of the original problem leveraging the BCU technique and random matrix theory is proposed.
First, P3 is transformed for better exposition based on the uplink-downlink duality [39] [45].
The instantaneous achievable rate in P3 is evaluated by the MIMO broadcast channel capacity.
The following Proposition 1 derives an implicit asymptotic expression of the objective function
in P3 such that the optimization is only dependent on statistical CSI which in this case is the
channel correlation matrices.
Proposition 1: In the large system regime, i.e., K → ∞ and K/Ns → β, the queue-weighted
downlink channel capacity in P3 is asymptotically equivalent to
D(Q, R1 , ..., RNt , Σu , Σb , p)


N
t
X
4
†
=
qn T log 1 + pn sn tr Σb BDFT
Rn BDFT Σb
n=1

1
M

n−1
X

†
pj sj Σb BDFT
Rj BDFT Σb

j=1

1 + en,j

!−1 
+I

 ,
(20)

13

where qi ’s are arranged in non-increasing order, and en,i is the unique solution of the following
equations.

†
Ri BDFT Σb
en,i = tr Σb BDFT

n−1
†
1 X pj sj Σb BDFT
Rj BDFT Σb
+I
M j=1
1 + en,j

!−1 
.

(21)

Proof: See Appendix B.
A. Convex Relaxation
Although the original MIP P3 is NP-hard, it can be relaxed to a multi-convex problem by
replacing the binary constraints with real-value constraints. The convex relaxation of an MIP is
a widely-used technique to achieve near-optimal solutions to the original problem [46] [47]. The
relaxed version of P3 is stated below.
P4:

maximize D(Q, R1 , ..., RNt , INt , Σb , w)
Σb ,w

s.t.,

Nt
X
n=1

wn ≤ P,

M
X

bi = K, 0 ≤ bi ≤ 1, ∀i,

(22)
(23)

i=1

where wn = pn sn , and D is defined in (20). Define the optimum solution of P4 as b∗i ’s and wn∗ ’s,
respectively. Then the scheduling decision is to schedule the beams and users corresponding to
the largest K b∗i ’s and Ns wn∗ ’s, respectively. It is observed that P4 is a multi-convex problem
[34] since the objective function is concave in both Σb and w. In view of this, the following
Algorithm 1, which bases upon the BCU technique is proposed.
The basic idea of the proposed BCU-based user and beam joint scheduling is that an iterative
method which cyclically optimizes user scheduling and beam selection with the other fixed
is guaranteed to converge to the global optimum of P4. In order to accelerate the iteration, an
iterative water filling approach which is based on [48] and deals with user scheduling is adopted.
Convergence of the proposed algorithm: The convergence to the global optimum is due to the
convergence results of the BCU algorithm [49]. The details of the proof is given in Appendix
C. It is found through simulations that the BCU-based algorithm normally converges after 23 iterations. Therefore, the computational complexity and convergence time are acceptable in
simulated scenarios.
Remark 2: The Algorithm 1 can solve the joint user scheduling and beam selection problem
based on statistical CSI. Therefore, it is applicable before channel estimations. After the system

14

Algorithm 1: BCU-Based Scheduling
(0)

1

Initialization: Σb = IM ;

2

Iteration: for t = 1 : T do

3

(0)

User scheduling update based on iterative water filling: ∀n ∈ [1, Nt ], ωn = P/Nt ;
for tw = 1 : Tw do

4

Compute for each n,
"
(t )
βn w

= tr

(t−1)
(t−1) †
Rn BDFT Σb
Σb BDFT



1
M

−1 #

(t −1)

wj w
Rj
j=1
1+en,j

Pn−1

+I

, where en,i is

the unique solution of the equations in (21);
5

Apply the classical water filling algorithm with water levels defined by β (tw )
γ

(tw )

=

arg max
P Nt

n=1

Nt
X

qn log 1 + γn βn(tw )



γn ≤P,γ≥0 n=1

Update ω as ω (tw ) = (1 − 1/M ) ω (tw −1) + (1/M )γ (tw )
6
7

8

(t )

(t −1)

if kωn w − ωn w

k <  then

w(t) = ω (tw ) , break;
(t)

Beam selection: Solve for Σb , which is the solution to the convex optimization
problem of P4 with w = w(t) .

9
10

11

(t)

(t−1)

Stopping criterion: if kw(t) − w(t−1) k < 1 and kΣb − Σb
wopt = ω (t) , Σb,opt =

(t)
Σb ,

k < 2 then

break;

Output: The scheduling user set is the users with the largest Nt values in wopt . The
selected beams are the ones with the largest K values in the diagonal entities of Σb,opt .

selects users and beams, the instantaneous channel estimations can be implemented and subsequently digital precoding and decoding can follow. This is in line with the multi-layer signal
processing concept proposed in, e.g., [4] [14], which proposes that the pre-beamforming should
be done based on channel statistics to save RF chains, complexity and system overhead.
V. I NCREMENTAL S ELECTION BASED M ETHOD FOR P3
Although the BCU-based algorithm is guaranteed to converge to the optimum of the relaxed
convex optimization problem, it still has high complexity due to the iterative algorithm design.

15

Therefore it may take a long time to converge. In this regard, an algorithm which selects users and
beams incrementally with low complexity is proposed. The key to the design of the algorithm is
to derive the incremental selection criterion. Thanks to the results in Proposition 1, the structure
of the asymptotic rates in (20) can be utilized to give such a criterion. Thereby, the incremental
selection algorithm is described in Algorithm 2, which assumes K = Ns 7 .
Algorithm 2: Incremental Greedy Scheduling (IGS)
1

Initialization: U = B = ∅; Uall = [1 : Nt ], Ball = [1 : M ]

2

Incremental Selection: for t = 1 : Ns do
Find the user-nt and the beam-bt that maximize:

3

(nt , bt ) =

arg max
n∈Uall \U,b∈Ball \B





qn log 1 +

P 
tr Σb R̄n Σb
Ns

P
1 X Ns Σb R̄j Σb
+I
M j∈U 1 + en,j

!−1 


†
where R̄j = BDFT
Rj BDFT ,

(Σb )i,i =


 1, for i ∈ B ∪ {b}


and


en,i = tr Σb R̄i Σb

0, else,

P
1 X Ns Σb R̄j Σb
+I
M j∈U 1 + en,j

!−1 
.

Update: U ∪ {nt } → U, B ∪ {bt } → B
Output: The scheduling user set is U, and the selected beam set is B.

4

Remark 3: Due to the successive interference cancellation (SIC) structure in the broadcast
channel queue-weighted capacity expression in Proposition 1, the rates of the users decoded
(selected in IGS) first will not be affected by the users decoded (selected in IGS) later. Therefore,
the IGS algorithm is viable because the earlier decisions are decoupled from later ones.
Complexity analysis: The IGS has a complexity of O(Ns Nt M ), because in each step it involves
an exhaustive search over O(M Nt ) possible user and beam combination, and there are Ns
7

The assumption is justified by arguing that maximum degree-of-freedom is achieved with K = Ns .

16

iterations. Compared with an exhaustive search over all user and beam subsets with complexity
 
Nt M
of N
, and the BCU-based scheduling which is difficult to quantize the complexity due to
K
s
the iterative design of the algorithm, the IGS has a relatively very low complexity.
Due to the greedy nature of the IGS algorithm, it may fail to find the optimum sets of users
and beams to schedule. Hence, it is better to have a worst-case performance bound of the IGS,
such that potentially arbitrarily bad solutions are excluded. Fortunately, this is the case for the
proposed IGS, due to the submodularity property [35] of the problem.
Theorem 2: Denote the queue weighted sum rate achieve by the IGS algorithm as DIGS , and
the global optimum as Dopt , then it is satisfied that
DIGS ≥ (1 − e−1 )Dopt .

(24)

Proof: The proof is based on the submodularity of the queue-weighted sum rate maximization problem in P3. Informally, the submodularity property indicates the problem has diminishing
returns, i.e., in this case the sum rate increase by scheduling a user or a beam is larger when
scheduled with a smaller user/beam set, i.e.,
D(U1 ∪ u) − D(U1 ) ≥ D(U2 ∪ u) − D(U2 ),

(25)

for any u ∈ Uall \(U1 ∪ U2 ) and U1 ⊆ U2 . This is easily validated since the same user will suffer
from more interference with a larger scheduled user set. There are two additional conditions
for submodularity, which is that the function should be nondecreasing and nonnegative. The
nondecreasing property, i.e.,
D(U ∪ u) ≥ D(U), ∀u ∈ Uall \U,

(26)

can be proved by arguing that at least, zero power can be allocated to the user or beam to obtain
equal performance without the user or beam since the objective function is a maximization over
all user and beam selection schemes. It should be noted that although the nonnegative condition
is easily validated, e.g., for the sum rate maximization or max-min rate maximization, it is not
met exactly for the PFS sum logarithm rate optimization. However, if we fix cn = 1, ∀n in (10),
the objective function is non-negative and thus the submodularity property is upheld.
Based on the submodularity, it can be proved that the IGS achieves a near-optimum performance. The remaining details of the proof is well-known in the literature and therefore omitted
for brevity [35], [50].

17

Remark 4: Although with Theorem 2, it is only proved that the IGS achieves at least about
60% throughput of the optimum scheme, the performance in reality is much better than that,
as which will be shown by simulations. Existing work which also utilizes the greedy algorithm
with submodularity property also agrees with this finding [50].
VI. S IMULATION R ESULTS
In this section, simulation results are presented. The channel model is as described in Section
II-B, where the number of MPCs for each user is Ln = 3, ∀n (including one LoS MPC), unless
stated otherwise. The ULA is used in the simulations and the amplitude of the LoS MPC is 10
times the one of the NLoS MPCs. The DoAs of the signals are generated from i.i.d. uniform
distributions. The antenna spacing d = λ/2, where λ denotes the carrier wavelength. In some of
the following cases where users’ pathlosses are not identical, the distances of users are generated
based on an i.i.d. uniform distributions from 30 to 200 meters and the pathloss γn is
 −γ
dn
,
γn =
d0

(27)

where γ = 2 which is in line with mm-wave channel measurements [51] and d0 is some reference
point distance. The regularized zero-forcing (RZF) precoder is adopted for system evaluation, i.e.,
−1
†
define Krzf = H̄ † H̄ + M αIM
. The RZF precoding matrix is expressed as Bd = ζKrzf H̄ ,
where ζ is a normalization scalar to fulfill the power constraint, and α is the regularization
factor [52]. Although RZF precoder is not the optimal coding scheme for Gaussian broadcast
channel (dirty-paper-coding with minimum-mean-square-error (MMSE) precoder is proved for
optimality but limited in reality due to high complexity), it can achieve full degree-of-freedom
(DoF) in the high signal-to-noise-ratio (SNR) region and is easy to implement. In the simulations,
α = Ns /ρ, where ρ is the receive SNR [52]. The user instantaneous rate is calculated by the
Shannon formula. The block fading model is adopted, where the channel stays constant for
10 time slots and evolves to another realization based on an i.i.d. distribution. The phase and
amplitude of each MPC is generated randomly. The simulation runs for 1000 such blocks and
calculate the time-averaged downlink transmission rates. The constants used in the Lyapunovdrift optimization are set to be V = Amax = 102 rmax , where rmax is the maximum rate of the
users. The , 1 and 2 in the stopping criterion in BCU-based algorithm are set to be 10−2 ρ/K.
In comparisons, the state-of-the-art interference aware beam selection scheme (IA-BS) in [21],

18

180
IA BS [21]
BCU-based
IGS
BDMA [29]

160

Sum rate (bits/s/Hz)

140
120
100
80
60
40
20
0
5

10

15

20

25

30

SNR (dB)

Fig. 2. Sum rate comparisons with identical user pathloss. The number of BS antennas is 256, the number of users is 100, the
number of scheduled users and beams are both 40.

and the BDMA scheme in [29] are also simulated. These schemes for comparisons are shown
to perform better than several other existing schemes, e.g., [22].
In Fig. 2, our proposed BCU-based algorithm and the IGS algorithm are compared with the
IA-BS algorithm [21] and the BDMA scheme [29]. Since the sum rates are considered, the
utility function in (9) is adopted. The scheduled user set in IA-BS is assumed to have the most
channel power. It is observed that by jointly considering user scheduling and beam selection, the
sum rate performance can be improved in the high-SNR regime. The reason is that in the highSNR regime, the interference is dominating the performance, and thus by jointly considering
the user scheduling and beam selection by the proposed schemes, the interference is better
suppressed. The BDMA scheme simply adopts a sequential approach which selects the beamuser pairs incrementally, and it builds on optimizing a sum-rate upper bound; both factors lead to
performance degradation. Thus, the resultant performance is not as good. Nonetheless, it should
be noted that the BDMA scheme is designed for multiple-antenna users and hence the interference
can be suppressed thereby whereas such effects are not captured in the presented simulations.
Therefore, we focus on the IA-BS scheme for comparisons in the following. Furthermore, the
figure also shows that the BCU-based and IGS algorithms achieve very similar performance in
this scenario.
Considering the utility functions with user-fairness considerations, e.g., the PFS utility function
in (10), the performance advantage of the proposed schemes is more obvious as shown in Fig. 3

19

0

100
IA BS [21]
BCU-based
IGS

-100

0
-50

Sum log rate

-200

Sum log rate

IA BS [21]
BCU-based
IGS

50

-300

-400

-100
-150
-200
-250

-500
-300
-INF

-350
5

10

15

20

25

30

5

10

SNR (dB)

15

20

25

30

SNR (dB)

Fig. 3. Sum log rate comparisons with (left) and without (right) non-identical user pathloss. The users with non-identical
pathloss are generated with distances i.i.d. uniformly distributed between 20 m to 200 m. The number of BS antennas is 256,
the number of users is 100, the number of scheduled users and beams are both 40.

(left). Based on the IA-BS, the beams with stronger channels are always preferred, corresponding
to users with small pathloss, resulting in ignorance of the fairness among users. In the proposed
joint scheduling schemes, the admission control in P2 utilizes a virtual queue to control the
fairness among users. Note that since there are always unscheduled users in the IA-BS due to
their small pathloss, the sum log rate of the IA-BS scheme in this case is negative infinity. In
Fig. 3 (right), the performance with identical pathloss is presented. Due to the fact that users
have identical large-scale fading, and hence equal probability to have good channels, the IA-BS
can achieve reasonably good fairness performance.
In order to show the performance loss of our proposed schemes compared with optimal
scheduling, an exhaustive search over all the feasible user and beam sets is conducted to solve
P3 and the optimal scheduling performance is obtained accordingly. Due to the prohibitive high
complexity of exhaustive search, we consider a small-scale problem where there are 8 users, 8
BS antennas and 4 scheduled users and beams. Nonetheless, it is found by many existing works,
e.g., [53] [54], that the impact of imperfect downlink scheduling decreases with the increase
of antenna dimension due to the channel hardening effect. Therefore, the relative performance
gap with a larger system dimension should be smaller, or at least similar with that in Fig. 4.
In Fig. 4, the left and right figures show sum-rate and sum-log-rate optimizations, respectively.
It is observed that in general the proposed schemes can achieve near-optimal performance. The
BCU-based scheme is shown to have better performance compared with the IGS, but with higher

20

35

12

10

25

Sum log rate

Sum rate (bits/s/Hz)

30

20
15
10

8

6

4

5
2
0
5

10

15

20

25

30

5

10

15

SNR (dB)

20

25

30

SNR (dB)

Fig. 4. Comparisons with optimal scheduling. The number of users is 8, the number of scheduled users and beams are both 4.
cn = 1, ∀n as in (10).

complexity. The performance bound which we prove for the IGS algorithm is also plotted in
the figure. The IGS scheme always performs better than the (1 − 1/e) bound. In the low SNR
regime, it is observed that the IA-BS scheme outperforms the IGS scheme, due to the reason that
the IA-BS scheme always selects the user and its corresponding beam which have the strongest
channel, and that in Fig. 4 we set cn = 1 in the log rate to ensure positive utilities and thus less
penalty on the unfairness among users is accounted for.
A throughput comparison with a stochastic traffic model is carried out and shown in Fig. 5
and 6. Instead of considering full-buffer greedy transmitting sources, each user’s arrival traffic
is modeled as a Bernoulli process, i.e., the arrival process for user-n
αn = bn rc , ∀n

(28)

where bn is i.i.d. Bernoulli distributed with expected mean values of pn , rc is a constant which
denotes approximately the service rate of each user, and rc =

Ns
Nt

log(1+η Nρs ) where η denotes the

approximate SNR loss coefficient introduced by interference. Therefore, based on this setting, the
mean values of bn can be regarded as the traffic intensity where pn = 0 denotes zero traffic, and
pn close to 1 denotes heavy traffic. The queuing dynamic is, with slightly abusing the notations,
qn (t + 1) = qn (t) − µ̃n (t) + αn (t),

(29)

where µ̃n (t) = min{qn (t), µn (t)} denotes the actual service rate taken into account of empty
queues. The throughput is calculated by averaging the sum actual service rate of each user µ̃n (t).

21

= 0.4

25

IA BS [21]
BCU-based
IGS

20

Throughput (bit/s/Hz)

20

Throughput (bit/s/Hz)

= 0.8

25

15

10

IA BS [21]
BCU-based
IGS

15

10

5

5

0

0
0

200

400

600

800

1000

Time slots

0

200

400

600

800

1000

Time slots

Fig. 5. Sample paths of the system average throughput evaluations. The number of BS antennas is 64, the number of users is
40, the number of scheduled users and beams are both 20. η = 0.4.

Our proposed Lyapunov-drift based schemes, i.e., BCU-based and IGS, can be easily adapted
to this traffic model, by replacing the virtual queues in (11) with (29) and eliminating the
admission control step in P2. We assume the IA-BS scheme [21] schedules Ns users with the
Ns -largest queue lengths and selects beams accordingly, which is in line with the methodology
that upper layers, e.g., medium-access-control (MAC) layer, schedules some users and push the
bits to physical layer. In comparisons, the proposed schemes jointly considers user scheduling
with traffic demands and beam selection. One sample path of the system is depicted in Fig. 5,
where the left and right sub-figures denote relatively low and heavy traffic, respectively. The
system average throughput is stationary after about 100 time slots, which equals about 100 ms
given the LTE numerologies where one time slot (transmission-time-interval) is one subframe
(1 ms). With this convergence time, it is found that the proposed schemes can effectively converge
to a reasonably good solution before the statistical CSI varies (usually at a time scale of several
seconds). In Fig. 6, the average throughput is compared among different scheduling schemes
under different traffic intensities. It is observed that the proposed schemes outperform the IA-BS
scheme when the system is with high traffic load. Note that when the system is not fully loaded,
i.e., traffic intensity is lower than about 0.6, the average throughput equals the sum arrival rate
and thus the performance advantage of the proposed schemes have not emerged. It is shown that
joint considerations of user scheduling and beam selection leads to better system throughput.

22

30
IA BS [21]
BCU-based
IGS

Throughput (bit/s/Hz)

25

20

15

10

0.4

0.5

0.6

0.7

0.8

0.9

1

Traffic intensity

Fig. 6. The number of BS antennas is 64, the number of users is 40, the number of scheduled users and beams are both 20.
η = 0.4.

VII. C ONCLUSIONS
In this paper, the BCU-based scheduling scheme and the IGS scheme are proposed to address
the join user scheduling and beam selection optimization problem in beam-based massive MIMO
systems based only on statistical CSI. The problem is formulated under the Lyapunov-drift
optimization framework. In order to solve the weighted rate maximization problem therein, the
proposed BCU-based scheme leverages the convex relaxation of the problem and adopts the BCU
technique with the iterative water-filling approach. It is proved that the BCU-based scheduling
scheme iteratively converges to the optimum of the relaxed problem. Due to its iterative algorithm
structure, relatively high complexity is required. Towards this end, the IGS algorithm is proposed
which is based on a greedy approach. Nevertheless, it is proved that the IGS scheme can achieve
performance within a multiplicative factor of (1−e−1 ) to the optimum. In simulations, it is shown
that the proposed schemes can achieve near-optimal performance and outperform the state-ofthe-art beam selection schemes, with utilities such as sum rate and proportional fairness. While
existing works focuses on the beam selection, which effectively strive to maximize the sum rate
performance, they are not optimized when the user scheduling and beam selection are jointly
considered especially when user fairness is taken into account. The performance bound we derive
for the IGS scheme is also shown to be well observed.

23

A PPENDIX A
P ROOF OF T HEOREM 1
We first briefly review the Lyapunov-drift approach, which is the main mathematical tool in
our proof. Define the Lyapunov function as
L(t) ,

1X 2
Q (t),
2 n n

(30)

and the Lyapunov drift as
∆(t) , E[L(t + 1) − L(t)|Q(t)].

Lemma 1: If there exist constants B and , which satisfy
X
Qn (t),
∆(t) ≤ B − 

(31)

(32)

n

then we have:
1) If  ≥ 0, then all queues are mean rate stable.
2) If  > 0, then

"
#
τ
X
1X
B
E
Qn (t) ≤ ,
lim sup

τ →∞ τ
t=1
n

(33)

and hence all queues are strongly stable.
Proof: The proof of Lemma 1 follows the standard procedure as in [32].
Given the queuing dynamics (11) and based on the definition in (31), we have
"
#

1X
∆(t) ≤ E
Q2n (t + 1) − Q2n (t) Q(t)
2 n
#
"

1X
=E
(Qn (t) − µ̃n (t) + an (t))2 − Q2n (t) Q(t)
2 n
"
#

X 1
1
=E
µ̃2n (t) − µ̃n (t)Qn (t) + a2n (t) + Qn (t)an (t) − an (t)µ̃n (t) Q(t)
2
2
n
"
#

X 1
1 2
2
≤E
µ (t) − µn (t)Qn (t) + an (t) + Qn (t)an (t) Q(t)
2 n
2
n
"
#
X µ2 (t) + a2 (t)
X
n
n
=E
Q(t) −
Qn (t)E [µn (t) − an (t)|Q(t)] .
(34)
2
n
n
Observing that
"
#
X µ2 (t) + a2 (t)

T2 
n
n
E
rn,max + A2max , B,
Q(t) ≤
2
2
n

(35)

24

where rn,max = log (1 + khn k2 P ) denotes the maximum transmission rate in one channel use
since rn,max is the channel capacity as if the user n was alone, it follows that
X
∆(t) ≤ B −
Qn (t)E [µn (t) − an (t)|Q(t)] .

(36)

n

Therefore, for any arrival rate inside the ergodic capacity region, since the scheduling problem
in P3 minimize the right-hand side of (36), the condition in Lemma 1 (32) is upheld with  ≥ 0,
i.e., all queues are mean rate stable. In order to show the throughput-optimality in (18), subtract
a term related to the utility function,
∆(t) − V E [U (a(t))| Q(t)] ≤ B −

X

Qn (t)E [µn (t)|Q(t)]

n

|

{z

}

Scheduling

"
+E

#
X

Qn (t)an (t) − V U (a(t)) Q(t) ,

(37)

n

|

{z

Admission control

}

It is observed that the admission control and scheduling problems in P2 and P3 are equivalent
to minimize the related terms in (37) as labeled. Therefore, given the solution to P2 and P3, the
left-hand side of (37) is less than the term on the right-hand side with any queue-independent
scheduling and admission control. Concretely,
∆(t) − V E [U (a(t))| Q(t)] ≤ B −

X

Qn (t)R̄n +

X

Qn (t)zn − V U (z),

(38)

n

n

where R̄n and zn denotes any queue-independent service rate and admission rate, respectively.
Taking expectations on both sides over Q(t), and taking the telescoping sum yields (assuming
Q(0) = 0 for better exposition),
τ −1

1X
E [Qn (t)] (R̄n − zn ) ≤ B + V
τ t=0

U

!
!
τ −1
1X
E [a(t)] − U (z)
τ t=0

(39)

Let z = R̄∗ which is the rate point in R that achieves the optimum utility function. Based on
the fact that all queues are mean rate stable as shown before, the left-hand side is non-negative,
it then follows that
∗

U (R̄ ) ≤ U

≤U

!
τ −1
1X
E [a(t)] + B/V
τ t=0
!
τ −1
1X
E [R(t)] + B/V
τ t=0

(40)

25

Let τ → ∞, take the lim sup and rearrange the terms yields the optimality condition in (18).
The inequality in (40) is based on the fact the queues are all mean rate stable the utility function
is non-decreasing. This completes the proof.
A PPENDIX B
P ROOF OF P ROPOSITION 1
Let qn = Qn (t) denote the virtual queue state at the scheduling time t, define
∆

G(Q, H, Σu , Σb , p) =

Nt
X

qn sn

n=1

=

=

Nt
X

T
X

qπ n sπ n

T
X

n=1

t=1

Nt
X

T
X

qπ n sπ n

Nt
X
n=1



P
†
det I + nj=1 Σb BDFT
hπj h†πj BDFT Σb pπj sπj


log
P
†
†
h
h
B
Σ
p
s
Σ
B
det I + n−1
b DFT πj πj DFT b πj πj
j=1


†
log 1 + h†πn BDFT Σb A−1 Σb BDFT
hπn pπn sπn

t=1

n=1

=

Rn (H(t), Σu , Σb , p)

t=1

qn

T
X



†
log 1 + h†n BDFT Σb A−1 Σb BDFT
hn pn sn

(41)

t=1

where
†
A = I + H̄[n−1] H̄[n−1]
,

H̄[n−1] =

√

†
pπj sπj Σb BDFT
[h1 , ..., hn−1 ] ,

(42)

and πi ∈ [1, ..., Nt ] is a permutation of the user index which satisfies qπ1 sπ1 ≥ ... ≥ qπNt sπNt
representing the decoding order in the dual uplink multiple-access-channel [40]. The last equality
is based on the fact that
x log det(I + Ax) = log det(I + Ax), ∀x ∈ {0, 1},

(43)

and without loss of generality, we assume qi ’s are arranged in non-increasing order. We invoke
[55, Theorem 1] which is stated at the end of the proof as Lemma 2 for reading convenience.

26

Denote the channel correlation matrix for user-n as Rn , and it yields

=
K→∞

−−−→

G(Q, H, Σu , Σb , p)
Nt
T

h
i
X
X
1
1
†
qn
log 1 + pn sn tr x†n Rn2 BDFT Σb A−1 Σb BDFT
Rn2 xn
n=1

t=1

Nt
X

T
X

n=1
K→∞

−−−→

Nt
X

qn



log 1 + pn sn tr

h

†
Rn BDFT Σb A−1
Σb BDFT

i

(44)

t=1





†
qn T log 1 + pn sn tr Σb BDFT
Rn BDFT Σb

n=1

n−1
†
1 X pj sj Σb BDFT
Rj BDFT Σb
+I
M j=1
1 + en,j

!−1 
 ,
(45)

where en,i is the unique solution of the following equations.

!−1 
n−1
†
X
1
pj sj Σb BDFT Rj BDFT Σb
†
.
en,i = tr Σb BDFT
Ri BDFT Σb
+I
M j=1
1 + en,j

(46)

The inequality (44) is based on [56, Lemma 14.2], and (45) is based on the following lemma.
Lemma 2: Let BN = XN† XN + SN with SN ∈ CN ×N Hermitian nonnegative definite and
1

XN ∈ Cn×N random. The ith column xi of XN† is xi = Ri2 yi , where the entries of yi ∈ Cri

are i.i.d. of zero mean, variance 1/N and have eighth-order moment of order O N14 . The
matrices Ri ’s are channel correlation matrices for each user, and QN ∈ CN ×N is deterministic.
Assume lim supN →∞ sup1≤i≤N kRi k < ∞ and let QN have uniformly bounded spectral norm
(with respect to N ). Define
mBN ,QN (z) =

1
trQN (BN − zIN )−1 .
N

(47)

Then, for z ∈ C\R+ , as n, N grow large with ratios βN,i = N/ri and β = N/n such that
0 < lim inf N βN ≤ lim supN βN < ∞ and 0 < lim inf N βN,i ≤ lim supN βN,i < ∞, we have that
mBN ,QN (z) − moBN ,QN (z) → 0

(48)

almost surely, with moBN ,QN (z) given by
moBN ,QN (z)

1
= trQN
N

n
1 X
Rj
+ SN − zIN
N j=1 1 + eN,j (z)

!−1
(49)

where the functions eN,j (z) form the unique solution of
1
eN,i (z) = trRi
N

n
1 X
Rj
+ SN − zIN
N j=1 1 + eN,j (z)

!−1
(50)

27

A PPENDIX C
P ROOF OF THE C ONVERGENCE OF THE BCU- BASED A LGORITHM
First, the proof of the iterative water filling approach in the user scheduling part is given.
Consider the user scheduling problem with beam selection fixed, i.e.,
P5:

maximize
w

Nt
X

qn log (1 + fn (w))

(51)

wn ≤ P,

(52)

n=1

s.t.,

Nt
X
n=1

Denote

†
fn (w) = fn (w1 , w2 , ..., wn ) = wn tr Σb BDFT
Rn BDFT Σb

n−1
1 X wj Rj
+I
M j=1 1 + en,j

!−1 
.

(53)

The key to the proof is to construct the equivalent optimization problem as stated below.
P6:

maximize

w(m), 0≤m≤Nt −1

Nt −1 X
Nt
1 X
qn log (1 + fn (w1 ([m + 1]Nt ), w2 ([m + 2]Nt ), ..., wn ([m + n]Nt )))
Nt m=0 n=1

(54)
s.t.,

Nt
X

wn (m) ≤ P, ∀m.

(55)

n=1

The reason that P5 and P6 are equivalent is straightforward due to the Shur-concavity of the
objective function of P6 [57]. Therefore, the solution of P6 is obtained at the point which satisfies
w(m) = w, ∀m.

(56)

Since P6 is concave in w(m), the BCU technique which cyclically optimizes w(m) with others
fixed is guaranteed to converge to the global optimum, which yields the same procedure as in
Algorithm 1 with some minor mathematical manipulations [58]. Therefore, we conclude that the
iterative water filling approach adopted in Algorithm 1 converges to the optimum in the user
scheduling step.
Next, it will be shown that the BCU technique which cyclically update user scheduling and
beam selection converges to the optimum. Based on [34], it is sufficient to check if the problem
satisfies the following two conditions:
•

The objective function, denoted by φ(x), is continuously differentiable in some neighborhood of every stationary point of φ(x).

28

•

For every k, 1 ≤ k ≤ n, φ(x) is a strictly concave function of xk , the other points xj ,
j 6= k, being arbitrarily chosen in their respective domains.

The above two conditions are easily met in this problem since D(Q, R1 , ..., RNt , INt , Σb , w)
is continuously differentiable in the whole domain and concave in Σb and w, respectively.
Therefore, the proposed BCU-based scheduling scheme is guaranteed to converge to the global
optimum.
R EFERENCES
[1] T. Marzetta, “Noncooperative cellular wireless with unlimited numbers of base station antennas,” IEEE Trans. Wireless
Commun., vol. 9, pp. 3590–3600, Nov 2010.
[2] F. Rusek, A. Lozano, and N. Jindal, “Mutual information of IID complex Gaussian signals on block Rayleigh-faded
channels,” IEEE Trans. Inform. Theory, vol. 58, pp. 331–340, Jan. 2012.
[3] E. G. Larsson, O. Edfors, F. Tufvesson, and T. L. Marzetta, “Massive MIMO for next generation wireless systems,” IEEE
Commun. Mag., vol. 52, pp. 186–195, Feb. 2014.
[4] R. W. Heath, N. Gonzlez-Prelcic, S. Rangan, W. Roh, and A. M. Sayeed, “An overview of signal processing techniques
for millimeter wave MIMO systems,” IEEE J. Sel. Top. Signal Process., vol. 10, pp. 436–453, Apr. 2016.
[5] S. Han, C. l. I, Z. Xu, and C. Rowell, “Large-scale antenna systems with hybrid analog and digital beamforming for
millimeter wave 5G,” IEEE Commun. Mag., vol. 53, pp. 186–194, Jan. 2015.
[6] Z. Jiang, S. Zhou, and Z. Niu, “Antenna-beam spatial transformation in C-RAN with large antenna arrays,” in IEEE Int.
Conf. Commun. (ICC) Workshops, May 2017.
[7] A. Alkhateeb, O. E. Ayach, G. Leus, and R. W. Heath, “Channel estimation and hybrid precoding for millimeter wave
cellular systems,” IEEE J. Sel. Top. Signal Process., vol. 8, pp. 831–846, Oct 2014.
[8] A. F. Molisch, V. V. Ratnam, S. Han, Z. Li, S. L. H. Nguyen, L. Li, and K. Haneda, “Hybrid beamforming for massive
MIMO: A survey,” IEEE Commun. Mag., vol. 55, pp. 134–141, Sep 2017.
[9] X. Gao, L. Dai, S. Han, C. L. I, and R. W. Heath, “Energy-efficient hybrid analog and digital precoding for mmwave
MIMO systems with large antenna arrays,” IEEE J. Select. Areas Commun., vol. 34, pp. 998–1009, Apr. 2016.
[10] J. Brady, N. Behdad, and A. M. Sayeed, “Beamspace MIMO for millimeter-wave communications: System architecture,
modeling, analysis, and measurements,” IEEE Trans. Antennas and Propagat., vol. 61, pp. 3814–3827, Jul. 2013.
[11] Z. Jiang, A. Molisch, G. Caire, and Z. Niu, “Achievable rates of FDD massive MIMO systems with spatial channel
correlation,” IEEE Trans. Wireless Commun., vol. 14, pp. 2868–2882, May 2015.
[12] Z. Jiang, S. Zhou, and Z. Niu, “On dimensionality loss in FDD massive MIMO systems,” in IEEE Conf. Wireless Commun.
Netw. (WCNC), pp. 399–404, Mar 2015.
[13] Z. Jiang, S. Zhou, R. Deng, Z. Niu, and S. Cao, “Pilot-data superposition for beam-based FDD massive MIMO downlinks,”
IEEE Commun. Letters, vol. 21, pp. 1357–1360, Jun 2017.
[14] A. Adhikary, J. Nam, J.-Y. Ahn, and G. Caire, “Joint spatial division and multiplexing: The large-scale array regime,”
IEEE Trans. Inform. Theory, vol. 59, pp. 6441–6463, Oct. 2013.
[15] A. F. Molisch and X. Zhang, “FFT-based hybrid antenna selection schemes for spatially correlated MIMO channels,” IEEE
Commun. Letters, vol. 8, pp. 36–38, Jan 2004.
[16] O. E. Ayach, S. Rajagopal, S. Abu-Surra, Z. Pi, and R. W. Heath, “Spatially sparse precoding in millimeter wave MIMO
systems,” IEEE Trans. Wireless Commun., vol. 13, pp. 1499–1513, Mar. 2014.

29

[17] T. Yoo and A. Goldsmith, “On the optimality of multiantenna broadcast scheduling using zero-forcing beamforming,”
IEEE J. Select. Areas Commun., vol. 24, pp. 528–541, Mar. 2006.
[18] H. Shirani-Mehr, G. Caire, and M. J. Neely, “MIMO downlink scheduling with non-perfect channel state knowledge,”
IEEE Trans. Commun., vol. 58, pp. 2055–2066, July 2010.
[19] H. Kim, K. Kim, Y. Han, and S. Yun, “A proportional fair scheduling for multicarrier transmission systems,” in IEEE Veh.
Tech. Conf. (VTC-Fall), vol. 1, pp. 409–413, Sep 2004.
[20] C. U. Bas, R. Wang, D. Psychoudakis, T. Henige, R. Monroe, J. Park, J. Zhang, and A. F. Molisch, “A real-time millimeterwave phased array MIMO channel sounder,” arXiv preprint arXiv:1703.05271, 2017.
[21] X. Gao, L. Dai, Z. Chen, Z. Wang, and Z. Zhang, “Near-optimal beam selection for beamspace mmwave massive MIMO
systems,” IEEE Commun. Letters, vol. 20, pp. 1054–1057, May 2016.
[22] P. V. Amadori and C. Masouros, “Low RF-complexity millimeter-wave beamspace-MIMO systems by beam selection,”
IEEE Trans. Commun., vol. 63, pp. 2212–2223, Jun 2015.
[23] Y.-C. Liang and F. P. S. Chin, “Downlink channel covariance matrix (DCCM) estimation and its applications in wireless
DS-CDMA systems,” IEEE J. Select. Areas Commun., vol. 19, pp. 222–232, Feb 2001.
[24] A. Adhikary, E. Al Safadi, M. Samimi, R. Wang, G. Caire, T. Rappaport, and A. Molisch, “Joint spatial division and
multiplexing for mm-wave channels,” IEEE J. Select. Areas Commun., vol. 32, pp. 1239–1255, Jun. 2014.
[25] A. Checko, H. L. Christiansen, Y. Yan, L. Scolari, G. Kardaras, M. S. Berger, and L. Dittmann, “Cloud RAN for mobile
networks – A technology overview,” IEEE Commun. Surveys Tuts., vol. 17, pp. 405–426, Firstquarter 2015.
[26] X. Gao, L. Dai, S. Han, C. L. I, and X. Wang, “Reliable beamspace channel estimation for millimeter-wave massive MIMO
systems with lens antenna array,” IEEE Trans. Wireless Commun., vol. 16, pp. 6010–6021, Sep 2017.
[27] Z. Jiang, S. Zhou, and Z. Niu, “A block coordinated update method for beam-based massive MIMO downlink scheduling
based on statistical CSI,” in Asia-Pacific Conference on Communications (APCC), Dec 2017.
[28] J. Wang, H. Zhu, L. Dai, N. J. Gomes, and J. Wang, “Low-complexity beam allocation for switched-beam based multiuser
massive MIMO systems,” IEEE Trans. Wireless Commun., vol. 15, pp. 8236–8248, Dec 2016.
[29] C. Sun, X. Gao, S. Jin, M. Matthaiou, Z. Ding, and C. Xiao, “Beam division multiple access transmission for massive
MIMO communications,” IEEE Trans. Commun., vol. 63, pp. 2170–2184, Jun 2015.
[30] L. You, X. Gao, G. Y. Li, X. G. Xia, and N. Ma, “BDMA for millimeter-wave/terahertz massive MIMO transmission with
per-beam synchronization,” IEEE J. Select. Areas Commun., vol. 35, pp. 1550–1563, Jul 2017.
[31] A. F. Molisch and M. Z. Win, “MIMO systems with antenna selection,” IEEE Microwave Mag., vol. 5, pp. 46–56, Mar
2004.
[32] M. J. Neely, “Stochastic network optimization with application to communication and queueing systems,” Synthesis Lectures
on Communication Networks, vol. 3, no. 1, pp. 1–211, 2010.
[33] A. Edelman and N. R. Rao, “Random matrix theory,” Acta Numerica, vol. 14, pp. 233–297, 2005.
[34] J. Warga, “Minimizing certain convex functions,” Journal of the Society for Industrial and Applied Mathematics, vol. 11,
no. 3, pp. 588–593, 1963.
[35] G. L. Nemhauser, L. A. Wolsey, and M. L. Fisher, “An analysis of approximations for maximizing submodular set
functions—I,” Mathematical Programming, vol. 14, no. 1, pp. 265–294, 1978.
[36] A. Molisch, “A generic model for MIMO wireless propagation channels in macro- and microcells,” IEEE Trans. Signal
Process., vol. 52, pp. 61–71, Jan. 2004.
[37] H. J. Kushner and P. A. Whiting, “Convergence of proportional-fair sharing algorithms under general conditions,” IEEE
Trans. Wireless Commun., vol. 3, pp. 1250–1259, Jul 2004.

30

[38] M. J. Neely, E. Modiano, and C. P. Li, “Fairness and optimal stochastic control for heterogeneous networks,” IEEE/ACM
Trans. Netw., vol. 16, pp. 396–409, Apr 2008.
[39] H. Weingarten, Y. Steinberg, and S. Shamai, “The capacity region of the Gaussian multiple-input multiple-output broadcast
channel,” IEEE Trans. Inform. Theory, vol. 52, pp. 3936–3964, Sep. 2006.
[40] W. Yu, “Uplink-downlink duality via minimax duality,” IEEE Trans. Inform. Theory, vol. 52, pp. 361–374, Feb. 2006.
[41] Z. Jiang, S. Zhou, and Z. Niu, “Capacity bounds of downlink network MIMO systems with inter-cluster interference,” in
IEEE Global Commun. Conf. (GLOBECOM), pp. 4612–4617, Dec. 2012.
[42] J. W. Lee, R. R. Mazumdar, and N. B. Shroff, “Downlink power allocation for multi-class CDMA wireless networks,” in
IEEE INFOCOM, vol. 3, pp. 1480–1489 vol.3, 2002.
[43] L. Georgiadis, M. J. Neely, and L. Tassiulas, Resource allocation and cross-layer control in wireless networks. Now
Publishers Inc, 2006.
[44] R. M. Karp, “Reducibility among combinatorial problems,” in Complexity of computer computations, pp. 85–103, Springer,
1972.
[45] W. Yu and J. Cioffi, “Sum capacity of Gaussian vector broadcast channels,” IEEE Trans. Inform. Theory, vol. 50, pp. 1875–
1892, Sep. 2004.
[46] S. E. Nai, W. Ser, Z. L. Yu, and H. Chen, “Beampattern synthesis for linear and planar arrays with antenna selection by
convex optimization,” IEEE Trans. Antennas and Propag., vol. 58, pp. 3923–3930, Dec 2010.
[47] A. Dua, K. Medepalli, and A. J. Paulraj, “Receive antenna selection in MIMO systems using convex optimization,” IEEE
Trans. Wireless Commun., vol. 5, pp. 2353–2357, Sep. 2006.
[48] N. Jindal, W. Rhee, S. Vishwanath, S. A. Jafar, and A. Goldsmith, “Sum power iterative water-filling for multi-antenna
Gaussian broadcast channels,” IEEE Trans. Inform. Theory, vol. 51, pp. 1570–1580, Apr 2005.
[49] P. Tseng, “Convergence of a block coordinate descent method for nondifferentiable minimization,” Journal of Optimization
Theory and Applications, vol. 109, no. 3, pp. 475–494, 2001.
[50] R. Kim, H. Lim, and B. Krishnamachari, “Prefetching-based data dissemination in vehicular cloud systems,” IEEE Trans.
Veh. Tech., vol. 65, pp. 292–306, Jan 2016.
[51] T. S. Rappaport, S. Sun, R. Mayzus, H. Zhao, Y. Azar, K. Wang, G. N. Wong, J. K. Schulz, M. Samimi, and F. Gutierrez,
“Millimeter wave mobile communications for 5G cellular: It will work!,” IEEE Access, vol. 1, pp. 335–349, 2013.
[52] C. Peel, B. Hochwald, and A. Swindlehurst, “A vector-perturbation technique for near-capacity multiantenna multiuser
communication-part I: Channel inversion and regularization,” IEEE Trans Commun., vol. 53, pp. 195–202, Jan. 2005.
[53] B. M. Hochwald, T. L. Marzetta, and V. Tarokh, “Multiple-antenna channel hardening and its implications for rate feedback
and scheduling,” IEEE Trans. Inform. Theory, vol. 50, pp. 1893–1909, Sep 2004.
[54] Z. Jiang, S. Zhou, and Z. Niu, “Dynamic channel acquisition in MU-MIMO,” IEEE Trans. Commun., vol. 62, pp. 4336–
4348, Dec. 2014.
[55] S. Wagner, R. Couillet, M. Debbah, and D. T. M. Slock, “Large system analysis of linear precoding in correlated MISO
broadcast channels under limited feedback,” IEEE Trans. Inform. Theory, vol. 58, pp. 4509–4537, Jul. 2012.
[56] R. Couillet and M. Debbah, Random matrix methods for wireless communications. Cambridge University Press, 2011.
[57] A. W. Marshall, I. Olkin, and B. C. Arnold, Inequalities: theory of majorization and its applications, vol. 143. Springer,
1979.
[58] M. Kobayashi and G. Caire, “An iterative water-filling algorithm for maximum weighted sum-rate of Gaussian MIMO-BC,”
IEEE J. Select. Areas Commun., vol. 24, pp. 1640–1646, Aug 2006.

