Stochastic Successive Convex Optimization for
Two-timescale Hybrid Precoding in Massive MIMO

arXiv:1801.08274v2 [] 29 Jan 2018

1

An Liu1 , Senior Member, IEEE, Vincent Lau1 , Fellow, IEEE and Min-Jian Zhao, Member, IEEE2
Department of Electronic and Computer Engineering, Hong Kong University of Science and Technology
2
College of Information Science and Electronic Engineering, Zhejiang University

Abstract—Hybrid precoding, which consists of an RF precoder
and a baseband precoder, is a popular precoding architecture
for massive MIMO due to its low hardware cost and power
consumption. In conventional hybrid precoding, both RF and
baseband precoders are adaptive to the real-time channel state
information (CSI). As a result, an individual RF precoder is
required for each subcarrier in wideband systems, leading to high
implementation cost. To overcome this issue, two-timescale hybrid
precoding (THP), which adapts the RF precoder to the channel
statistics, has been proposed. Since the channel statistics are
approximately the same over different subcarriers, only a single
RF precoder is required in THP. Despite the advantages of THP,
there lacks a unified and efficient algorithm for its optimization
due to the non-convex and stochastic nature of the problem.
Based on stochastic successive convex approximation (SSCA), we
propose an online algorithmic framework called SSCA-THP for
general THP optimization problems, in which the hybrid precoder is updated by solving a quadratic surrogate optimization
problem whenever a new channel sample is obtained. Then we
prove the convergence of SSCA-THP to stationary points. Finally,
we apply SSCA-THP to solve three important THP optimization
problems and verify its advantages over existing solutions.
Index Terms—Massive MIMO, Two-timescale Hybrid Precoding, Successive Convex Approximation

I. I NTRODUCTION
Massive MIMO can significantly improve the spectrum
efficiency of wireless systems. The conventional fully digital precoding requires one RF chain for each antenna, and
thus induces huge hardware cost and power consumption for
massive MIMO. As a result, hybrid precoding, where a highdimensional RF precoder is connected to a low-dimensional
baseband precoder with a limited number of RF chains,
has been proposed to reduce the hardware cost and power
consumption of massive MIMO base station (BS).
The early works on this topic focus on studying fasttimescale hybrid precoding (FHP), where both RF and baseband precoders are adaptive to the real-time channel state
information (CSI). For example, in [1], a sparse precoding and
combining algorithm based on orthogonal matching pursuit
is proposed for single-user mmWave systems. In [2], a lowcomplexity FHP scheme for multiuser massive MIMO systems
is proposed. A limited feedback hybrid precoding scheme is
also proposed in [3] for multi-user mmWave systems. One
disadvantage of FHP is that the number of RF precoders
has to increase with the number of subcarriers in wideband
systems such as orthogonal frequency-division multiple access
(OFDMA) system (since the CSI is different on different
subcarriers), leading to high implementation cost [4], [5].
Moreover, as the real-time full CSI is required at the BS, the
CSI signaling overhead is large.

To overcome the above disadvantages of FHP, a twotimescale hybrid precoding (THP) scheme is proposed in [6],
[7]. In THP, the RF precoder is adaptive to the channel
statistics1 to achieve the array gain, and the baseband precoder
is adaptive to the low-dimensional effective channel to achieve
the spatial multiplexing gain. THP has several advantages.
Since the channel statistics are approximately the same on
different subcarriers [8], THP only needs one RF precoder to
cover all subcarriers, which significantly reduces the hardware
cost. Moreover, THP reduces the CSI signaling overhead
because it does not require knowledge of the real-time highdimensional CSI. Therefore, THP can achieve a better tradeoff
between the implementation cost and performance, making it
more attractive in practice [4], [5].
The optimal THP design depends on the RF precoding
structure and the specific application scenario. There are two
major RF precoding structures: the fully-connected structure
where each antenna is connected to all the RF chains, and
the partially-connected structure where each antenna is only
connected to a single RF chain [9]. For each structure, there
are two methods to implement the dynamic RF precoder.
In the codebook-based method, the RF precoder is chosen
from a pre-determined codebook [1], [3], [6], while in the
dynamic-phase-shifter-based (DPS-based) method, the phase
of each element of the RF precoder can be quantized and
adjusted individually [2], [5]. Under different RF precoding
structures/implementations, the constraint on the RF precoder
is different and thus the optimal THP design is also different.
Moreover, in different application scenarios, the optimization
objectives can also be quite different. For example, for besteffort services, we may want to maximize the throughput or
proportional fairness (PFS) utility under a total power constraint. For applications with a fixed throughput requirement,
such as video streaming, we may want to minimize the transmit power subject to an individual throughput requirement.
Therefore, it is important to develop a systematic solution
framework to optimize the THP design for a wide range of
applications (i.e., with a general objective function) under
different RF precoding structure/implementation constraints.
Unfortunately, the optimization of THP is quite challenging due to the non-convex stochastic optimization problem
involved. The existing solutions are usually heuristic and
only suitable for one application scenario under a specific
RF precoding structure constraint. For example, the minimum weighted throughput maximization (MWTM) problem
1 In this paper, channel statistics refers to the moments or distribution of
the channel fading realizations.

2

is approximately solved in [6] under the fully-connected and
DFT-based RF precoding structure (i.e., the RF precoding
codebook forms a DFT matrix). Specifically, the average data
rate is first replaced with a closed-form lower bound based
on knowledge of the channel covariance matrices, and then
the resulting approximate problem is solved by semidefinite
relaxation (SDR). However, the lower bound may become
loose when different user clusters have overlapped angle of
departure (AoD) intervals or the SNR is low. In [4], the
average sum-rate maximization problem is solved using the
sample average approximation (SAA) method, again under the
fully-connected and codebook-based RF precoding structure.
However, it is known that SAA has high complexity and is
not suitable for online implementation. The Signal-to-leakageand-noise ratio (SLNR) maximization problem is considered
in [5] under the fully-connected and DPS-based RF precoding
structure. But the algorithm in [5] does not consider the
fairness issue.
Note that all the above methods are offline methods, which
require a channel sample collection phase to collect a large
number of channel samples (to estimate the channel covariance
matrices or construct the SAA functions) before computing
the optimized hybrid precoder, and thus its performance at the
channel sample collection phase is limited. In [10], a bestresponse-based (BRB) algorithm is proposed for solving general stochastic non-convex multi-agent optimization problems.
The BRB algorithm is an online algorithm based on stochastic
successive convex optimization (SSCA). However, it only
works when the objective function contains expectations but
the constraint can be represented by a deterministic convex set.
In many application scenarios, such as the MWTM problem
considered in [6], there are stochastic non-convex constraints
(i.e., the constraint functions are also non-convex and involve
expectations over the random states) involved, which are
difficult to deal with.
In this paper, we propose an online algorithmic framework
called SSCA-THP to solve a general THP optimization problem without explicit knowledge of channel statistics. The main
contributions are summarized as follows.
• A general THP optimization formulation: We propose
a general THP optimization formulation with a general
smooth objective function, which can be applied to various application scenarios under different RF precoding
structure/implementation constraints.
• An online algorithmic framework based on SSCA and
its convergence proof: We propose an online algorithmic framework called SSCA-THP to solve the general
THP optimization problem with stochastic non-convex
constraints, and establish its convergence to stationary
points. At each iteration of SSCA-THP, quadratic surrogate functions are constructed for both objective and
constraint functions based on a new channel realization
and the current iterate. Then the next iterate is updated
by solving the resulting quadratic optimization problem
using a low-complexity Lagrange dual method. SSCATHP has several advantages over existing algorithms.
First, it is an online algorithm, meaning that the RF
precoder is updated whenever a (potentially outdated)

DPS
FHP
MWTM
MM
SSCA
SAA
THP

Dynamic phase shifter
Fast-timescale hybrid precoding
Minimum weighted throughput maximization
Majorization-minimization
Stochastic successive convex approximation
Sample average approximation
Two-timescale hybrid precoding

Table I: List of abbreviations.

Figure 1: Massive MIMO downlink with THP

channel sample is obtained. As a result, it can achieve
a better overall performance than the offline counterpart.
Second, it only requires outdated full CSI samples and
thus is more robust to signaling latency in practical wireless networks. Third, the quadratic optimization problem
at each iteration can be efficiently solved by the Lagrange
dual method, which has very low complexity. Finally,
SSCA-THP provides a systematic solution for the design
of THP, which opens the door to solving the more difficult
THP optimization problems that occur in practice.
• Specific SSCA-THP algorithm design for important
applications: We apply SSCA-THP to solve several
important THP optimization problems in massive MIMO.
Simulations verify the advantages of the proposed algorithmic framework over existing baseline solutions.
The rest of the paper is organized as follows. In Section
II, we present the system model for the massive MIMO
downlink with THP, various implementation methods for the
RF precoder, and the general THP optimization formulation.
The SSCA-THP algorithm and the convergence analysis are
presented in Section III and IV, respectively. Section V applies SSCA-THP to solve several important THP optimization
problems. Finally, the conclusion is given in Section VI.
Notations: |S| denotes the cardinality of a set S. Diag (a)
represents a diagonal matrix whose diagonal elements form
the vector a. For a matrix M , Diag (M ) denotes a
vector consisting of the diagonal elements of M and
[M ]i,j denotes the (i, j)-th element of M . Let M =
BlockDiag (M 1 , M 2 , ..., M n ) denote a block diagonal matrix
with the i-th block given by M i , and S [M ] , M + M H .
Let ◦ denote the Hadamard product and R[M ] denote the real
part of a complex matrix M .
II. S YSTEM M ODEL AND P ROBLEM F ORMULATION
A. Multi-user Massive MIMO Downlink with THP
Consider a multi-user massive MIMO downlink system with
one BS serving K single-antenna users, as illustrated in Fig.
1. For clarity, we focus on a narrowband system with flat
block fading channel, but the proposed algorithm can be easily

3

Figure 3: Fully-connected and partially-connected structures in RF
precoder

Figure 2: Timeline (frame structure) of SSCA-THP algorithm

modified to cover the wideband system as well. The BS is
equipped with M antennas and S transmit RF chains, where
K ≤ S  M . Hybrid precoding is employed to support
simultaneous transmissions to the K users with limited RF
chains at the BS. In this case, the transmit signal vector
for user k is given by F g k sk , where F ∈ CM ×S is the
RF precoder, and g k ∈ CS×1 and sk ∼ CN (0, 1) are the
baseband precoding vector and the data symbol for user k,
respectively. The RF precoder F is usually implemented using
an RF phase shifting network [11]. Hence, all
√ elements of
F have equal magnitude, i.e., [F ]i,j = √1M e −1θi,j , where
θi,j is the phase of the (i, j)-th element of F . Under hybrid
precoding, the received signal for user k is given by
X√
√
H
yk = pk hH
pi F g i si + zk ,
(1)
k F g k sk + hk
i6=k

where hk ∈ CM is the channel of user k, pk is the transmit
power allocated to user k, and zk ∼ CN (0, 1) is the additive
white Gaussian noise (AWGN).
In this paper, we consider a THP scheme, whose time line
is illustrated in Fig. 2. The time domain is divided into superframes. Each super-frame is further divided into L frames, and
each frame consists of Tf time slots. The channel statistics
(distribution) are assumed to be constant within each superframe, and the channel state H is assumed to be constant
within each time slot. We assume that the BS can obtain
the real-time effective CSI HF at each time slot, and one
(possibly outdated) channel sample H at each frame. In our
design, the BS is not required to have explicit knowledge of
the channel statistics. By observing one channel sample at
each frame, the proposed algorithm can automatically learn
the channel statistics (in an implicit way). Specifically, the
RF precoder F is only updated once per frame based on the
(possibly outdated) channel sample to achieve massive MIMO
array gain with reduced implementation cost. On the other
hand, the baseband precoder G = [g 1 , ..., g K ] is adaptive to
the real-time effective CSI HF ∈ CK×S to achieve the spatial
H
multiplexing gain, where H = [h1 , ..., hK ] ∈ CK×M is the
composite downlink channel. We consider a regularized zeroforcing (RZF) baseband precoder [12]:

−1
G = F H H H HF F H H H + αI
Λ1/2 ,
(2)
whereh α is the regularization
factor, Λ
=
i
−1
−1
is used to normalize the
Diag kg 1 k , ..., kg K k

precoding vectors
 F g k ’s, and g k is−1the k-th column of
H
H
HF F H H H + αI
.
G , FF H
Although the baseband precoder G is adaptive to the
instantaneous effective CSI, as in (2), we assume that the regT
ularization factor α and the power allocation p = [p1 , ..., pK ]
are adaptive to the channel statistics only. This is because
in the massive MIMO regime, the system tends to behave
like a deterministic system and thus the gain of adapting the
power allocation and regularization factor according to the
instantaneous CSI is small [13], [14]. Similar assumptions
have also been made in [5], [6], [15] to achieve a good
compromise between performance and complexity.
B. Various Implementation Methods for RF Precoder
Various implementation methods for the RF precoder F
have been proposed in the literature to achieve different
tradeoffs between the performance, complexity and power
consumption. Basically, there are two major RF precoding
structures: the fully-connected structure and the partiallyconnected structure. For each structure, there are two common
methods to dynamically adjust the RF precoder. In this subsection, we shall discuss the existing typical implementations
for the RF precoder and their pros and cons.
One major challenge for the optimization of the RF precoder
is that it has a discrete implementation constraint; e.g., the
phase shifter cannot take a continuous value in practice and
has to be quantized into discrete values. In this subsection,
we will also discuss two techniques to convert the discrete RF
precoder into continuous variables to make the optimization
of the RF precoder tractable.
1) Fully-connected RF Precoding Structure: In this case,
each RF chain is connected to every antenna through phase
shifters and RF adders. To be more specific, if a BS has M
antennas and S RF chains, M S phase shifters and M RF
adders will be required to implement the RF precoder, as
illustrated in Fig. 3-(a). There are two methods to dynamically
adjust the RF precoder, as elaborated below.
DPS-based RF Precoder: In this method, each phase
shifter is quantized using B bits. Theneach phase shifter θ
i,j
2π (2B −1)
2π
can take a value from the discrete set 0, 2B , ...,
.
2B
Since the phase of each element of the RF precoder can
be quantized and adjusted individually, the DPS-based RF
precoder can achieve a good performance. However, the requirement of a large number of dynamic phase shifters and RF
adders increases the hardware cost and power consumption.

4

The DPS-based RF precoder can be represented by a phase
vector θ ∈ RM S whose ((j − 1) M + i)-th element is θi,j ,
and the optimization of the DPS-based RF precoder is a
discrete optimization problem due to the discrete constraint on
θi,j . To make the problem tractable, we first ignore the discrete
constraint and treat θi,j as a continuous variable. After finding
the optimized phase vector θ ∗ , we project it onto the feasible
set to obtain the final phase vector θ̂:

∗
θ̂i,j =
argmin
θi,j
mod 2π − θ , (3)
θ∈{0,...,2π(2B −1)/2B }

∀i, j, where mod denotes the modulo operation. With only
B = 3 bits, the performance loss due to the quantization effect
is already small [16].
Codebook-based RF Precoder: In this case, each column
of the RF precoder can only be selected from a finite-size
codebook F = {c1 , ..., cN } with |F| = N code vectors.
The RF precoding codebook F is usually chosen to be the
array response vectors at the BS. For example, when a uniform
linear array (ULA) is used, the array response vectors form a
DFT matrix. The codebook-based RF precoder can be easily
implemented using a static precoder at the RF domain (using a
static phase shifting network) together with an RF switch [6].
However, the performance of the codebook-based RF precoder
is in general worse than that of the DPS-based RF precoder.
The codebook-based RF precoder can be represented by a
selection matrix as F = CD, where C = [c1 , ..., cN ] ∈
CM ×N is the codebook matrix, P
and D , Diag (d) with
T
N
N
d = [d1 , ..., dN ] ∈ {0, 1} and i=1 di = S is a selection
matrix. The optimization of the codebook-based RF precoder
is an integer optimization problem due to the constraint
di ∈
PN
{0, 1} , ∀i = 1, ..., N . Note that the constraints i=1 di = S
and di ∈ {0, 1} , ∀i are equivalent to the following constraints:
N
X

di = S, di ∈ [0, 1] , kdk0 ≤ S.

(4)

i=1

To make the problem tractable, we approximate the l0 -norm
kdk0 using a smooth function as [17]




N
X
1
di
N
/ log 1 +
, d ∈ [0, 1] , (5)
kdk0 ≈
log 1 +


i=1
where the smooth parameter  > 0 can be used to control the
approximation error. A smaller  leads to smaller approximation error but a less smooth function. Then we can replace the
constraint kdk0 ≤ S with a non-convex sparse constraint




N
X
1
di
log 1 +
/ log 1 +
≤ S,
(6)


i=1
which can be handled using the SSCA-THP algorithm. Note
that the sparse constraint (6) has been widely used in compressive sensing to impose sparsity on sparse signals [17]. After
finding the optimized selection vector d∗ , we project it onto
the feasible set to obtain the final solution d̂:
N

d̂ = argmin kd − d∗ k , s.t. d ∈ {0, 1} , kdk0 = S.
d

Specifically, d̂ is a binary vector whose S non-zero elements
are located at the S largest elements of d∗ . Thanks to the

sparse constraint on d∗ , d∗ is usually close to a binary vector
with S non-zero elements, and thus the performance loss due
to the above projection is very small.
2) Partially-connected RF Precoding Structure: In this
case, each RF chain is connected to a sub-array of antennas
via phase shifters and each antenna is only connected to
a single RF chain [9], as illustrated in Fig. 3-(b). Such a
partially-connected structure presents a block diagonal form
as F = BlockDiag (f 1 , f 2 , ..., f S ), where f s ∈ CM/S is
the s-th precoding vector corresponding to the s-th RF chain.
Note that we have assumed that M is divisible by S for easy
illustration. The partially-connected structure only requires
a total number of M phase shifters (instead of M S phase
shifters in the fully-connected case). As a result, it has much
lower power consumption and hardware complexity compared
to the fully-connected case. However, the performance of the
partially-connected RF precoder is compromised.
The DPS-based or codebook-based method can be used to
adjust the partially-connected RF precoder. In the DPS-based
method, the partially-connected RF precoder can be represented by a phase vector θ ∈ RM whose ((j − 1) M/S + i)th element is θi,j for i = (j − 1) M/S + 1, ..., jM/S;
j = 1, ..., S. Similarly, to make the problem tractable, we
first ignore the discrete constraint on θi,j and then project the
resulting solution onto the discrete set using (3).
In the codebook-based method, the s-th RF precoding
vector f s is selected from a finite-size codebook Fs =
{c1 , ..., cN } with |Fs | = N code vectors. Let C s =
M
[cs,1 , ..., cs,N ] ∈ C S ×N denote the codebook matrix for
f s . Then the s-th RF precoding vector f s can be represented as f s = C s D s , where D s P, Diag (ds ) with
T
N
N
ds = [ds,1 , ..., ds,N ] ∈ {0, 1} and
i=1 ds,i = 1 is a
selection matrix. For convenience, let d = [d1 , ..., dS ] denote
the overall selection vector. Similarly, we can use the “sparse”
technique to relax the integer constraints ds,i ∈ {0, 1} , ∀s, i
PN
to non-convex smooth
constraints
d
= 1, ds,i ∈
i=1


 s,i
PN
ds,i
1
[0, 1] , i=1 log 1 +  / log 1 +  ≤ 1.
C. Achievable Data Rate
Under different implementation methods, the RF precoder
F is represented by different parameters with different dimensions. For convenience, we use φ as a unified notation
to denote the RF precoding parameter. For example, in the
DPS-based method, φ = θ and its dimension is M S and
M for the fully-connected and partially-connected structures,
respectively. In the codebook-based method, φ = d and its
dimension is N and N S for the fully-connected and partiallyconnected structures, respectively.
For given RF precoding parameter φ, power allocation p,
the RZF regularization factor α in (2) and channel realization
H, the instantaneous achievable data rate of user k is


2
H
p
h
F
g
k
k
k


rk (φ, p, α; H) = log 1 + P
.
2
H
p
h
F
g
+
1
i
i
k
i6=k
Note that F is a function of φ and G is a function of φ, p, α
and H. Therefore, we explicitly express rk as a function of

5

φ, p, α which depends on the random channel state H. The
average data rate of user k is
rk (φ, p, α) = E [rk (φ, p, α; H)] .
For
convenience,
define
r (φ, p, α)
,
[r1 (φ, p, α) , ..., rK (φ, p, α)]T as the average data rate
vector.
D. THP Optimization Formulation
Consider the following general optimization formulation for
the design of THP:
min f0 (x) , h0 (r (φ, p, α) , x)

(7)

x∈X

s.t.fi (x) , hi (r (φ, p, α) , x) ≤ 0, i = 1, ...., m,
iT
h
where x = φT , pT , α, β T
is called the THP variable,
β ∈ Rnβ , with dimension nβ , is an additional optimization
variable that does not affect the average data rate vector r, and
X ⊆ Rn is a convex set. The motivation for introducing an
additional optimization variable β is to provide extra flexibility
in the formulation to cover more application scenarios, as
will be shown in Example 4 at the end of this subsection.
Note that both the dimension n and the set X depend on the
implementation method for the RF precoder. For example, for
the codebook-based RF precoder, we have
n
o
N
K
X = x : d ∈ [0, 1] ; p ∈ [0, p̃] ; α > α̃; β ∈ B , (8)

where ε > 0 is a small number used to avoid the singularity
at rk = 0 [19].
PKThis is an instance of Problem
PK (7) with
h0 (r, x) = − k=1 log (ε + rk ), h1 (r, x) = k=1 pk − P
and B = ∅.
Example 3 (Power minimization with individual QoS requirements [20]). In this example, the THP variable x is designed
to minimize the average transmit power subject to individual
QoS requirements as follows:
min

x∈X

K
X

pk , s.t. rk (φ, p, α) ≥ γk , ∀k,

(11)

k=1

where each user has an individual QoS requirement in terms
of the average data rate constraint rk (φ, p, α) ≥ γk , and the
constant γk ≥ 0 is the target rate for
user k. This is an instance
PK
of Problem (7) with h0 (r, x) = k=1 pk , hk (r, x) = γk −
rk , ∀k and B = ∅.
Example 4 (MWTM [6]). In this example, the THP variable
x is designed to maximize the minimum (weighted) average
data rate of users as follows:
K
X
1
pk ≤ P,
max min
rk (φ, p, α) , s.t.
x∈X
k wk

(12)

k=1

where wk > 0 is the weight for user k, which can be used
to provide a differential QoS for different users. (12) is not
an instance of Problem (7) because the objective function is
non-smooth. However, by introducing an auxiliary variable β,
we can convert (12) to an instance of Problem (7) as:

where p̃ > 0 is used to ensure that the transmit power
is bounded, α̃ > 0 is a small number to ensure that the
calculation of the matrix inverse in (2) is always numerically
stable, B = {β : βi ∈ Bi , i = 1, ..., nβ }, and Bi is a convex
region in R. Note that, without loss of generality, all the
coupled constraints on x, such as the sparse constraint (6) on
the codebook-based RF precoding parameter d, are included in
the explicit constraints fi (x) ≤ 0, i = 1, ...., m. As a result, X
has a decoupled form: X = {x : xi ∈ Xi , i = 1, ..., n}, where
Xi is a convex region in R. We assume that the functions
hi (r, x) , i = 0, ..., m are continuously differentiable (and
possibly non-convex) functions of (r, x).
Problem (7) embraces many applications as special cases.
In the following, we give some important examples.

There are several challenges to solve Problem (7). First,
the average data rates rk (φ, p, α)’s are neither convex nor
concave, and have no closed-form expressions. Moreover, the
presence of stochastic non-convex constraints further complicates Problem (7). In the next section, we shall propose an
efficient algorithm based on the SSCA method, called SSCATHP, to find a stationary point of Problem (7).

Example 1 (Sum throughput maximization [4]). The sum
throughput maximization problem is formulated as:

III. S TOCHASTIC S UCCESSIVE C ONVEX A PPROXIMATION
FOR THP O PTIMIZATION

max
x∈X

K
X
k=1

rk (φ, p, α) , s.t.

K
X

(9)

k=1

Example 2 (PFS [18]). PFS is a widely used utility function in
wireless resource optimization. The PFS utility maximization
problem is formulated as:

x∈X

K
X
k=1

log (ε + rk (φ, p, α)) , s.t.

x∈X

pk ≤ P,

(13)

k=1

with h0 (r, x) P = −β, hk (r, x) = wk β − rk , ∀k,
K
hK+1 (r, x) = k=1 pk − P and B = {β : β ≥ 0}.

A. The SSCA-THP Algorithm
pk ≤ P,

where P is the total power constraint at the BS.PThis is
K
an instance of Problem (7) with h0 (r, x) = − k=1 rk ,
PK
h1 (r, x) =
k=1 pk − P and B = ∅ (i.e., there is no
additional variable β).

max

max β, s.t. rk (φ, p, α) ≥ wk β, ∀k;

K
X

K
X
k=1

pk ≤ P,

(10)

At each iteration, the THP variable x is updated by solving
a quadratic optimization problem obtained by replacing the
objective and constraint functions fi (x) , i = 0, ..., m with
their quadratic surrogate functions f¯il (x) , i = 0, ..., m.
Algorithm 1 summarizes the key steps of the proposed
SSCA-THP algorithm. Specifically, at iteration l, a new realization of the random channel state H l is obtained in Step
1 and the surrogate functions f¯il (x) , ∀i are updated based on
H l and the current iterate xl as


T

2
f¯il (x) = hi r̂ l , xl + uli
x − xl + τi x − xl ,
(14)

6

Algorithm 1 SSCA-THP Algorithm


Input: γ l , ρl .
Initialize: x0 ∈ X ; u−1
= 0, ∀i, l = 0.
i
Step 1: Obtain a channel sample H l within frame l.
Update the surrogate functions f¯il (x) , ∀i using (14).
Step 2: Solve (17) to obtain the optimal solution ν ◦ , x◦ .
If ν ◦ ≤ 0 (Problem (16) is feasible)
Solve (16) to obtain x̄l . //Objective update
Else
Let x̄l = x◦ . //Feasible update
End if
Step 3: Update xl+1 according to (18).
Step 4: Let l = l + 1 and return to Step 1.

B. Efficient Solutions for Quadratic Optimization Subproblems



l T
where τi > 0 is a constant; r̂ l = r̂1l , ..., r̂K
, with
P
l
l
j
l
l
r
φ
,
p
,
α
;
H
/l,
is
the
sample
average
r̂kl =
j=1 k


approximations for rk φl , pl , αl ; uli is an approximation for

the gradient ∇fi xl , which is updated recursively as

uli = 1 − ρl ul−1
+ ρl ûli ,
i
with u−1 = 0, where ρl ∈ (0, 1] is a sequence to be properly
chosen and






ûli = Jr xl ; H l ∇r hi r̂ l , xl + ∇x hi r̂ l , xl , (15)


where Jr xl ; H l
is the Jacobian matrix of
the instantaneous rate vector r (φ, p, α; H)
,
[r1 (φ, p, α; H) , ..., rK (φ, p, α; H)]T and its expression is
derived in Appendix A, ∇r hi and ∇x hi are the gradients of
hi w.r.t. the average rate vector r and the THP variable x,
respectively. The surrogate function f¯il (x) can be viewed as
a convex approximation of fi (x) in a local domain around
xl .
In Step 2, the optimal solution x̄l of the following problem
is solved:
x̄l = argmin f¯0l (x)

In this subsection, we propose efficient solutions for the
quadratic optimization subproblems in (16) and (17) based
on the Lagrange dual method. The reasons for using the
Lagrange dual method are as follows. First, for given Lagrange
multipliers (which are also called dual variables), the problem
of minimizing the Lagrange function has a unique and closedform solution. Second, the number of primal variables x is
usually much larger than the number of constraints (dual
variables) in the massive MIMO regime. Therefore, the optimal dual variables can be solved much more efficiently than
directly solving the optimal primal variables.
In the following, we show how to use the Lagrange dual
method to solve the subproblem in (16). Since both subproblems have the same form (i.e., both are strictly convex
and quadratic optimization problems), the solution for (17) is
similar and is omitted for conciseness.
The Lagrange function for (16) is
L (x, λ) = f¯0l (x) +
l

λi f¯il (x)

i=1

=

n
X

a (λ) x2i + bi (λ) xi + c (λ) , x ∈ X ,

i=1
T

where λ = [λ1 , ..., λm ] ∈ Rm
+ are the Lagrange multipliers,
a (λ) =
bi (λ) =

m
X
j=0
m
X

λj τj ,

λj ulj,i − 2τj xli ,

j=0

c (λ) =

m
X

 

T
λj hj r̂ l , xl − ulj xl + τj xl



,

λ0 = 1, and xli and ulj,i are the i-th element of xl and ulj ,
respectively. The dual function for (16) is

(16)

g l (λ) = min Ll (x, λ) .
x∈X

s.t. f¯il (x) ≤ 0, i = 1, ...., m,

2

j=0

x∈X

(19)

And the corresponding dual problem is

which is a convex approximation of (7). Note that Problem
(16) is not necessarily feasible. If Problem (16) turns out to
be infeasible, the optimal solution x̄l of the following convex
problem is solved:
x̄l = argmin ν

m
X

max g l (λ) .
λ≥0

(20)

The minimization problem in (19) can be decomposed into
n independent subproblems as
min a (λ) x2i + bi (λ) xi , i = 1, ..., n,

(17)

x∈X ,ν

s.t. f¯il (x) ≤ ν, i = 1, ...., m,
which minimizes the constraint functions.
Given x̄l in one of the above two cases, x is updated in
Step 3 according to

xl+1 = 1 − γ l xl + γ l x̄l ,
(18)
where γ l ∈ (0, 1] is a sequence to be properly chosen. Then the
above iteration (Steps 1 to 3) is carried out until convergence.

xi ∈Xi

which have the following closed-form solutions:


bi (λ)
◦
, ∀i,
xi (λ) = PXi −
2a (λ)

(21)

where PXi [·] denotes the one-dimensional projection on
to the convex set Xi . On the other hand, the dual
function g (λ) is concave and it can be verified that
 l ◦
T
l
f¯1 (x (λ)) , , ..., f¯m
(x◦ (λ)) is a subgradient of g (λ) at
λ. Hence, the standard subgradient-based methods such as
the subgradient algorithm in [21] or the ellipsoid method in

7

[22] can be used to solve the optimal solution λ◦ of the dual
problem in (20). Then the optimal primal solution of (16) is
given by x◦ (λ◦ ).

4) Let x∗F be any stationary point of the following feasibility problem:
min ν

(22)

x∈X ,ν

C. Implementation Consideration
At the beginning of each super-frame, the BS resets the
SSCA-THP algorithm with an initial THP variable x0 . Then
the THP variable x is updated once every frame. Therefore,
each frame corresponds to an iteration
al  in the SSCA-THP
 T T
T

T
gorithm. Specifically, let xl = φl , pl , αl , β l
denote the THP variable used during the l-th frame. At time
slot t in the l-th frame, the BS first acquires the effective
channel H (t) F l , where H (t) is the channel state at time
slot t, and F l is the RF precoder corresponding to φl . Then
it calculates the baseband precoder G (t) from H (t) F l and
pl , αl using (2). At the end of the l-th frame, the BS obtains a
channel sample H l and updates the THP variable x by solving
a simple quadratic optimization problem, where the updated
THP variable xl+1 will be used in the (l + 1)-th frame. Then
the same procedure is carried out in the next frame.
Remark 1. The proposed SSCA-THP algorithm exploits some
unique properties of hybrid beamforming massive MIMO systems to improve the performance and reduce the complexity.
For example, the optimization variables and constraints in
Section II-B are specifically designed for hybrid beamforming
with different RF precoding structures. By imposing the sparse
constraint in (6), the property of limited RF chains in massive
MIMO is also exploited to improve the performance over
the existing semidefinite relaxation (SDR) method in [4] for
the codebook-based RF precoder. Finally, the structure of the
average data rate function rk (θ, p, α) w.r.t. the phase vector
θ is exploited to design quadratic surrogate functions which
enables low-complexity and fast-convergent algorithm design.
IV. C ONVERGENCE A NALYSIS
In this section, we establish the local convergence of SSCATHP to a stationary point. There are several challenges in
the convergence proof for SSCA-THP. First, we need to show
that at every limiting point, all constraints are satisfied, which
is non-trivial since SSCA-THP may oscillate between the
feasible update and objective update. Moreover, the limiting
point is obtained by averaging over all the previous outputs
from either feasible updates or objective updates, which makes
it difficult to show that the limiting point is a stationary
point of the original Problem (7). In this subsection, we will
overcome these challenges and establish the convergence of
SSCA-THP. To prove the convergence of SSCA-THP, we need
to make the following assumptions on the problem structure.
Assumption 1 (Assumptions on Problem (7)).
1) hi (r, x) , i = 0, ..., m are continuously differentiable
functions of (r, x).
2) For any x ∈ X , the functions hi (r (φ, p, α) , x) , i =
0, ..., m, their derivative, and their second-order derivative w.r.t. r and x are uniformly bounded.
3) H l , l = 0, 1, ... are uniformly bounded w.p.1.

s.t. fi (x) ≤ ν, ∀i = 1, ...., m.
We assume that fi (x∗F ) ≤ 0, i = 1, ..., m.
The first assumption is standard and is satisfied for a large
class of problems. In practice, the channel sample is always
bounded, and thus the second assumption is satisfied. The third
assumption ensures that Problem (7) is feasible. If there is a
stationary point x∗F which is not feasible, then Algorithm 1
may get stuck at this stationary point x∗F . Therefore, the third
assumption is necessary for the algorithm to converge to a
feasible point of the problem.
Besides Assumption 1, the sequence of parameters
{ρt } , {γ t } needs to satisfy the following conditions.
Assumption 2 (Assumptions on {ρt } , {γ t }).
2
P
P
1) ρl → 0, l ρl = ∞, l ρl < ∞, liml→∞ ρl l−1/2 <
∞.
2
P
P
2) γ l → 0, l γ l = ∞, l γ l < ∞,
3) liml→∞ γ l /ρl = 0.
With Assumptions 1 and 2, we can prove two key lemmas
that will eventually lead to the final convergence result. The
first lemma proves the convergence of surrogate functions.
Lemma 1 (Convergence of surrogate functions). Suppose
Assumptions
1 and 2 are satisfied. Consider a subsequence
 l ∞
x j j=1 converging to a limiting point x∗ , and define functions
fˆi (x) , hi (r (φ∗ , p∗ , α∗ ) , x∗ )
2

+ ∇fi (x∗ ) (x − x∗ ) + τi kx − x∗ k , ∀i,
which satisfy fˆi (x∗ ) = fi (x∗ ) and ∇fˆi (x∗ ) = ∇fi (x∗ ) , ∀i.
Then, almost surely, we have
l
lim f¯i j (x) = fˆi (x) , ∀x ∈ X .

j→∞

(23)

Please refer to Appendix B for the proof. To state the
convergence result, we need to introduce the concept of Slater
condition for the converged surrogate functions.
Slater condition for
 the∞converged surrogate functions:
Given a subsequence xlj j=1 converging to a limiting point
x∗ and letting fˆi (x) , ∀i be the converged surrogate functions
as defined in Lemma 1, we say that the Slater condition is
satisfied at x∗ if there exists x ∈ intX such that
fˆi (x) < 0, ∀i = 1, ..., m.
A similar Slater condition is also assumed in [23] to prove
the convergence of a deterministic majorization-minimization
(MM) algorithm with non-convex constraints.
Before the introduction of the main convergence theorem,
we give the second key lemma.

8

l→∞

Average Tx power (dB)


∞
Lemma 2. Let xl l=1 denote the sequence of iterates
generated by Algorithm 1. We have

lim sup max fi xl ≤ 0, w.p.1.
i∈{1,...,m}

lim x̄l − xl = 0, w.p.1.

The lemma states that the algorithm will converge to the
feasible region, and the gap between x̄l and xl converges to
zero, almost surely. Please refer to Appendix C for the proof.
Theorem 1 (Convergence of Algorithm 1). SupposeAssump∞
tions 1 and 2 are satisfied. For any subsequence xlj j=1
∗
converging to a limit point x , if the Slater condition is
satisfied at x∗ , then x∗ is a stationary point of Problem (7)
almost surely.
Please refer to Appendix D for the proof.
V. A PPLICATIONS AND N UMERICAL VALIDATION
In this section, we shall apply the proposed SSCA-THP to
solve the first three example problems described in Section
II. As in [5], we adopt a geometry-based channel model with
a half-wavelength space ULA for simulations.
PNp The channel
vector of user k can be expressed as hk = i=1
αk,i a (ϕk,i ),
where a (ϕ) is the array response vector, ϕk,i ’s are Laplacian distributed
with an angle spread σAS = 10, αk,i ∼

2
2
CN 0, σk,i
, σk,i
are randomly generated from an exponenPNp 2
tial distribution and normalized such that
i=1 σk,i = gk ,
and gk represents the path gain of user k. Unless otherwise
specified, we consider M = 64 antennas, S = 16 RF chains
and Np = 6 channel paths. The path gains gk ’s are uniformly
generated between -10 dB and 10 dB. The DPS-based RF
precoder is considered in Example 1 and 2, while both the
DPS-based and codebook-based RF precoders are considered
in Example 3. We compare the performance of the SSCA-THP
with the following baseline algorithms.
Baseline 1 (SAA): This is the sample average approximation algorithm. Specifically, after applying SAA, the problem
becomes a deterministic non-convex optimization problem,
which is then solved using the deterministic successive convex
approximation method [23].
Baseline 2 (SLNR-max): This is the SLNR maximization
algorithm in [5].
Baseline 3 (JSDM): This is the the joint spatial division and
multiplexing scheme in [15].
Baseline 4 (CB): This is the CB algorithm in [6], which
combines the deterministic approximation and bisection methods.
Both SSCA-THP and SAA can be used to solve a general
THP optimization problem, while the SLNR-max/JSDM is
more suitable for the sum throughput maximization in Example 1, and the CB algorithm can be used to solve the power
minimization problem in Example 3. Both the SLNR-max and
CB algorithms only work for the fully-connected RF precoding
structure. For fair comparison, we focus on the fully-connected
structure in the simulations. All baseline algorithms belong to
the offline method, which requires a channel sample collection
phase to construct the SAA functions (baseline 1) or estimate

Maximum constraint (bps/Hz)

l→∞

8
7
6
5
4

0

50

100
Number of iterations

150

200

0

50

100
Number of iterations

150

200

0.6
0.4
0.2
0
−0.2

Figure 4: Convergence of the SSCA-THP

the channel covariance matrices (baseline 2 and 3). We assume
that one super-frame has L = 1000 frames and the first 200
frames serve as the channel sample collection phase for the
baseline algorithms. For fair comparison, the proposed SSCATHP is also terminated after 200 frames (iterations). The
performance is obtained by averaging over the last 800 frames
of the super-frame. Note that if we considered the overall performance averaged over the entire super-frame, the proposed
SSCA-THP would achieve an even larger performance gain
over the baseline algorithms, which perform poorly during the
channel sample collection phase.
A. Convergence of the Proposed SSCA-THP
We use Example 3 to illustrate the fast convergence of
SSCA-THP. Specifically, there are K = 8 users and the target
average rate for all users is set to be the same as γk = 2
bps/Hz. Consider the DFT-based RF precoder. In Fig. 4, we
plot the objective function (average transmit power) and the
maximum constraint function (target average rate minus the
minimum achieved average rate of users) versus the iteration
number, respectively. It can be seen that SSCA-THP quickly
converges to a point with all target average rates satisfied with
high accuracy.
B. Sum Throughput Maximization
In Figs. 5 and 6, we plot the sum throughput versus the
number of users K and the number of channel paths Np ,
respectively. The transmit power is set to be P = 10 dB.
For comparison, we also plot the sum throughput of the
fully digital RZF beamforming [12] and the FHP algorithm
based on alternating optimization (AO-FHP) in [16]. The
proposed SSCA-THP achieves better performance than the
existing THP algorithms (SLNR-max and JSDM). Moreover,
as the number of users increases, the performance gap between
SSCA-THP and SLNR-max/JSDM increases. Although SAA
achieves similar performance to SSCA-THP after the channel
sample collection phase, it has poor performance during the
channel sample collection phase. Moreover, the per iteration
complexity of SAA is much higher than that of SSCA-THP

9

80

26

Sum spectrum efficiency (bps/Hz)

75
70
65
60
55

22

Fully digital

20

Proposed, 3 bits
Proposed & SAA,
unquantized

50
45
SLNR−max, 3 bits

SLNR-max,
unquantized

16
14

SLNR-max,
2, 3, 4 bits

6

8

10
12
Number of users

Proposed, unquantized
SAA, unquantized
SLNR-max, unquantized
Proposed, quantized
SLNR-max, quantized

10

SLNR−max & JSDM,
unquantized

35

8

14

6

16

6

7

8

9

10

11

12

13

14

15

16

Number of users

Figure 5: Sum throughput (bps/Hz) versus the number of users K

Figure 7: PFS utility versus the number of users K

58

30

56

52
Proposed & SAA,
unquantized
Proposed, 3 bits

48
46

Proposed
SAA
SLNR−max
Fully Digital
AO−FHP
JSDM
Proposed, 3 bits
SLNR−max, 3bits

Total transmit power (dB)

54

50

Proposed, DFT codebook
SAA, DFT codebook
CB, DFT codebook
Proposed, unquantized DPS
Proposed, quantized DPS

25

AO−FHP

Fully digital
Sum spectrum efficiency (bps/Hz)

18

12

40

30

Proposed,
2, 3, 4 bits

Proposed & SAA,
unquantized

24

AO−FHP

PFS utility

Proposed
SAA
SLNR−max
Fully Digital
AO−FHP
JSDM
Proposed, 3 bits
SLNR−max, 3bits

20

Prop. & SAA,
DFT
15

10
Proposed, DPS
unquantized

5

Proposed, DPS
1, 2, 3, 4 bits

44
SLNR−max & JSDM,
unquantized

42

CB, DFT

0
SLNR−max, 3 bits

40

2

4

6
8
Number of channel paths

10

12

−5

1

2

3
4
5
Target average data rate per user (bps/Hz)

6

Figure 6: Sum throughput versus the number of channel paths Np

Figure 8: Average transmit power versus the target average rate

(CPU time: 10.00 s versus 0.06 s). The performance gap
between the THP with statistical RF precoder and the FHP
with real-time RF precoder (or fully digital RZF) is smaller
when the number of users/channel paths is smaller. This is
consistent with the analysis in [4]. Finally, it can be seen that
with only three-bit quantization, the performance is already
very close to the case without quantization.

K = 8 users. For any given target average rate γ, the proposed
SSCA-THP with the DPS-based RF precoder achieves the
lowest transmit power. The performance of the DPS-based RF
precoder degrades as the number of quantization bits B for
each phase decreases. When B = 1, the performance of the
DPS-based RF precoder is similar to that of the DFT-based RF
precoder optimized using the SSCA-THP or SAA algorithms,
which is still much better than the DFT-based RF precoder
optimized using the CB algorithm.

C. Proportional Fairness
In Fig. 7, we plot the PFS utility versus the number of
users K. The transmit power is set to be P = 10 dB. Similar
results to those in Fig. 5 can be observed. Moreover, when
considering the PFS utility, the performance gap between
SSCA-THP and SLNR-max is much larger since the fairness
issue is not considered in the SLNR-max algorithm. Note that
the PFS utility of SLNR-max decreases with the number of
users. This is because, without considering the fairness, the
minimum throughput of the users becomes much smaller as
the number of users increases.
D. Power Minimization
In Fig. 8, we plot the average transmit power versus the
target average rate requirement γk = γ, ∀k for a system with

E. Complexity Comparison
The complexity of SSCA-THP is dominated

 by the call
l
culation of the Jacobian matrix Jr x ; H , which has
complexity order O (M KS), and the quadratic optimization
subproblems
 in (16) and (17), which has complexity order
O M m2 , as explained below. For given Lagrange multipliers, the complexity order of calculating the closed-form primal
solution in (21) is O (M ). Using the ellipsoid method, the
number of iterations required to achieve a given convergence
accuracy  for the dual problem (20) is O m2 log (1/)
[22]. Hence, the per-iteration
complexity of SSCA-THP is

O M m2 log 1 + M KS . In Table II, we compare the complexity order of SSCA-THP with SLNR-max in [5] and JSDM

10

SSCA-THP

O LI M log

SLNR-max
1


+ M KS



O LI M 2 S + M 3



JSDM 
3

O M

Table II: Comparison of the complexity order for different algorithms.

in [15] for the sum throughput maximization problem (i.e.,
m = 1), where LI is the total number of iterations for SSCATHP and Algorithm 1 in [5], respectively. The complexity
order of both SLNR-max and JSDM increases with the number
of BS antennas M according to M 3 because they involve
singular value decomposition (SVD) for the M × M channel
covariance matrix. On the other hand, the complexity order
of SSCA-THP only increases linearly with M , thanks to the
closed-form solution in (21) for fixed Lagrange multipliers.

VI. C ONCLUSION
In this paper, we first propose a general optimization formulation (7) for the design of THP in massive
MIMO, which is applicable to different RF precoding structures/implementations and a wide range of application scenarios. Then we propose an online algorithmic framework called
SSCA-THP to solve this general THP optimization problem.
Specifically, at each iteration, quadratic surrogate functions are
constructed for both objective and constraint functions based
on a new channel sample. Then the next iterate is updated
by solving the resulting quadratic optimization problem. We
prove the convergence of SSCA-THP to stationary points.
To the best of our knowledge, SSCA-THP is the first online
and provably convergent algorithm to handle the general nonconvex stochastic constraints considered in (7). Finally, we
apply SSCA-THP to solve three important THP optimization
problems and verify its advantages.

A PPENDIX
A. Jacobian Matrix of Instantaneous Rate
1) Jacobian Matrix for the Fully-connected DPS-based RF
MS
Precoder:
In this
and x =
h
iT case, we have φ = θ = R
T
T
T
θ , p , α, β
. We first define some useful notations:
i
h
H
H
Ak,i = 2H H
F S HF F hk hk GP i Λ H F F
h
i
− 2S hk hH
GP
ΛH
F,
i
F
k
h
i


H
E i = 2H H
F S HF F GI i H F F − 2S GI i H F F ,
h
i
H
H
2
ek,i = H F F F H hk hH
,
k F F H F P iΛ
i,i



H

H

−1

where H F = BH with B = HF F H + αI
, and
P i (I i ) denotes a K × K matrix with [P i ]i,i = pi ([I i ]i,i =
1) and all other elements being zero. Then using the matrix

calculus, it can be shown that the gradients of rk (θ, p, α; H)
w.r.t. θ, p and α are respectively given by
P
P θ
θ
i6=k ak,i
i ak,i
∇θ rk (θ, p, α; H) =
−
(24)
Γ
Γ
P −k p
P kp
i6=k ak,i
i ak,i
−
,
(25)
∇p rk (θ, p, α; H) =
Γ
Γ
P −k α
P kα
i6=k ak,i
i ak,i
∇α rk (θ, p, α; H) =
−
,
(26)
Γk
Γ−k
where
P
i6=k

Γk

P

=

pi hH
k F gi

i

pi hH
k F gi

2

+ 1,

Γ−k

=

2

+ 1,

√

√

aθk,i = Vec R −1F ∗ ◦ Ak,i − ek,i Vec R −1F ∗ ◦ E i ,
 H

apk,i = Diag G hk hH
k GΛi ,
h 
i
h
i
H
H
aα
− 2R hH
k,i = ek,i 2R T r GI i BG
k GP i Λi BG hk .

Therefore, for given channel state H, the Jacobian matrix of
the instantaneous rate vector r (θ, p, α; H) w.r.t. x is


∇θ r1 ∇θ r2 · · · ∇θ rK
 ∇p r1 ∇p r2 · · · ∇p rK 

(27)
Jr (x; H) = 
 ∇α r1 ∇α r2 · · · ∇α rK  ,
0
0
0
0
where the bottom submatrix is zero because ∇β rk = 0, ∀k.
Note that we have omitted (θ, p, α; H) in the gradient expressions for simplicity of notation.
2) Jacobian Matrix for the Fully-connected CodebookN
based RF Precoder: In this case, we have φ = d = [0, 1]
h
iT
and x = dT , pT , α, β T . Using the matrix calculus, it can
be shown that the gradients of rk (d, p, α; H) w.r.t. p and α
are given by (25) and (26), respectively; and the gradient of
rk (d, p, α; H) w.r.t. d is
P
P d
d
i6=k ak,i
i ak,i
∇d rk (θ, p, α; H) =
−
,
Γ
Γk




e
H
where adk,i = − 21 Diag C H Ak,i C + k,i
2 Diag C E i C .
Finally, for given channel state H, the Jacobian matrix of the
instantaneous rate vector r (d, p, α; H) w.r.t. x is given by
(27) with ∇θ rk , ∀k replaced by ∇d rk , ∀k.
The Jacobian matrix of the instantaneous rate vector for the
partially-connected structure can be obtained similarly. The
details are omitted for conciseness.
B. Proof of Lemma 1
The proof relies on the following lemma.
Lemma 3. Under Assumption 2, we have


lim f¯il xl − fi xl = 0,
l→∞


lim ∇f¯il xl − ∇fi xl = 0,
l→∞


lim f¯il1 xl1 − f¯il2 xl2 ≤ C xl1 − xl2 .
l1 ,l2 →∞

for i = 0, ..., m w.p.1., where C > 0 is some constant.

(28)
(29)
(30)

11

Proof: It follows from the law of large numbers and the
central limit theorem that
 
1
l a.s. l
l
l
r̂ → r , E r̂ − r = O √ ,
(31)
l


where r l = r φl , pl , αl . Then (28) follows from (31).
On the other hand, (29) is a consequence of [24], Lemma
1. It is easy to verify that the technical conditions (a), (b), (d)
and (e) therein are satisfied. In the following, we prove that
l
condition (c) in [24], Lemma 1 is also satisfied.
Let ∇ξ hi =


ˆ l hi = ∇ξ hi r̂ l , xl for ξ ∈ {r, x}. Let
∇ξ hi r l , xl and ∇
h  ξ i

Jr xl = E Jr xl , H l denote the Jacobian matrix of the


average rate vector r φl , pl , αl at point xl . Then we have

 l
l
∇fi xl = Jr xl ∇r hi + ∇x hi .

Case 2: x◦ is the optimal solution of the following convex
optimization problem:
min ν
s.t. fˆi (x) ≤ ν, i = 1, ...., m.
Since the Slater condition is satisfied (by choosing a sufficiently large ν, we can always find a point x ∈ X such that
fˆi (x) < ν, i = 1, ...., m), the KKT condition of the problem
(35) implies that there exist λ1 , ..., λm such that
X
λi ∇fˆi (x◦ ) = 0,
i

1−

(32)
λi




 

ˆ lr hi − ∇lr hi
E ûli − ∇fi xl ≤ E Jr xl , H l ∇
ˆ lx hi − ∇lx hi
+E ∇
 


1
a
(33)
=O r̂ l − r l = O √ ,
l
where (33-a) holds because ∇hi are Lipschitz continuous and Jr xl , H l are bounded w.p.1. From (33) and

 
P∞ l −0.5
P∞
< ∞, we have l=1 ρl E ûli − ∇fi xl <
l=1 ρ l
∞, which implies that the technical condition (c) in [24],
Lemma 1 is satisfied.
Finally, (30) follows from the Lipschitz continuity of hi .
This completes the proof.
n From oLemma 3 and
n (31),o the families of functions
l
f¯i j (x) converge to fˆi (x) almost surely.
C. Proof of Lemma 2

1. We first prove lim supl→∞ f xl ≤ 0 w.p.1., where
f (x) = maxi∈{1,...,m}fi (x).
Let T = l : f xl ≥  for any  > 0. We show that T
is a finite set by contradiction.
Suppose T is infinite. We first show that
lim inf l∈T ,l→∞ x̄l − xl > 0 by contradiction. Suppose
lim inf l∈T ,l→∞ x̄l − xl
= 0. Then there exists a
subsequence lj ∈ T such that limj→∞ x̄lj − xlj = 0. Let
x◦ denote a limiting point of the subsequence xlj , and let
fˆi (x) , ∀i be the converged surrogate functions as defined in
Lemma 1. According to the update rule of Algorithm 1, there
are two cases.
Case 1: x◦ is the optimal solution of the following convex
optimization problem:

x

X

λi = 0,

i

It follows from (32) and (31) that

min fˆ0 (x)

(35)

x,ν

(34)

s.t. fˆi (x) ≤ 0, i = 1, ...., m.
In this case, we have f (x◦ ) = maxi∈{1,...,m} fˆi (x◦ ) ≤ 0,
which contradicts the definition of T .



fˆi (x◦ ) ≤ ν, ∀i = 1, ..., m,

fˆi (x◦ ) − ν = 0, ∀i = 1, ..., m.

(36)

It follows from Lemma 1 and (36) that x◦ also satisfies
the KKT condition of Problem (22). By Assumption 1, we
have fi (x◦ ) ≤ 0, i = 1, ..., m, which again contradicts the
definition of T .
Therefore, lim inf l∈T ,l→∞ x̄l − xl > 0; i.e., there exists
a sufficiently large l such that
0

0

x̄l − xl ≥  , ∀l ∈ T ,
0

(37)

0

where  > 0 is some constant and T = T ∩ {l ≥ l }.
Define function f¯l (x) = maxi∈{1,...,m} f¯il (x). From the
definition of f¯il (x) in (14), f¯il (x) is strongly convex, and
thus
2



∇T f¯il xl dl ≤ −η dl + f¯il x̄l − f¯il xl ,
(38)
where dl = x̄l − xl , and η > 0 is some constant. From
Assumption 1, the gradient of fi (x) is Lipschitz continuous,
and thus there exists Lf > 0 such that



2 l 2
fi xl+1 ≤ fi xl + γ l ∇T fi xl dl + Lf γ l
d
2




2
= f xl + Lf γ l
dl + fi xl − f xl



+ γ l ∇T f¯il xl + ∇T fi xl − ∇T f¯il xl dl
2
a



≤ f xl + fi xl − f xl − ηγ l dl



+ γ l f¯il x̄l − f¯il xl + o γ l
2


≤ f xl − ηγ l dl + o γ l , ∀i = 1, ..., m,
(39)


where o γ l means that liml→∞ o γ l /γ l = 0. In (39-a),
we used (38) and liml→∞ ∇T fi xl − ∇T f¯il xl  = 0,
and the last inequality follows from fi xl
≤
f xl , lim infl→∞ f xl − f¯il x̄l
≥
0, and
liml→∞ fi xl − f¯il xl
= 0. Since (39) holds for
all i = 1, ..., m, by choosing a sufficiently large l , we have
2


f xl+1 − f xl ≤ −γ l η dl
0

0

≤ −γ l η , ∀l ∈ T ,

(40)

12

for some η > 0. Moreover, from Assumption 1, the directional
derivative of f (x) is uniformly bounded, and thus there exists
a constant C such that


0
(41)
f xl+1 − f xl ≤ C xl+1 − xl ≤ C γ l ,
0

for some C > 0. Finally, it follows from (40) and (41) that

f xl ≤ 2, ∀l ≥ l .
(42)
Since (42) is true
for any  > 0, it follows that

lim supl→∞ f xl ≤ 0.
2. Then we prove that liml→∞ x̄l − xl = 0, w.p.1.
2.1: We first prove that lim inf l→∞ x̄l − xl = 0 w.p.1.
Note that the feasible problem in (17) is strictly convex, and
thus the solution is uniquely given by x̄l . Therefore, when
 a
feasible update is performed at iteration l, we have f¯l x̄l ≥ 0
and
x̄l = argmin f¯0l (x)
x


s.t. f¯il (x) ≤ f¯l x̄l , i = 1, ...., m.
As a result, x̄l can be expressed in a unified way as
x̄l = argmin f¯0l (x)
x

s.t. f¯il (x) ≤ ν l , i = 1, ...., m,
where ν = 0 when an objective update is performed and
ν l = f¯l x̄l when a feasible
update is performed.
Since


liml→∞ f¯l xl − f xl = 0, f¯l x̄l ≤ f¯l xl , and we
have proved that lim supl→∞ f xl ≤ 0, it follows that
liml→∞ ν l = 0. Let x̂l denote the projection of xl on
to the feasible set of Problem (43).
 Then it follows from

liml→∞ ν l = 0, lim supl→∞ f¯l xl =
lim supt→∞ f xl ≤

0, and the strong convexity of f¯l xl that
(44)

From the definition of f¯0l (x) in (14), f¯0l (x) is uniformly
strongly convex, and thus
2



∇T f¯0l xl dl ≤ −η dl + f¯0l x̄l − f¯0l xl
 
2

l
= −η dl + f¯0l x̄l − f¯0 x̂l
 

l
+ f¯0 x̂l − f¯0l xl
≤ −η dl

2

+ e (l) ,

It follows from (46) that
l
X

2

f0 xl − f0 xl0 ≤ −ηχ2
γj ,
j=l0

2
P∞
which, in view of  j=l0 γ j
= ∞, contradicts
the boundedness of f0 xl . Therefore it must be
lim inf l→∞ x̄l − xl = 0 w.p.1.
2.2: Then we prove that lim supl→∞ x̄l − xl = 0 w.p.1.
We first prove a useful lemma.
x̄l1 − x̄l2 ≤ L̂ xl1 − xl2 + e (l1 , l2 ) ,
where liml1 ,l2 →∞ e (l1 , l2 ) = 0.

l

lim xl − x̂l = 0.

probability. Then we can find a realization such that dl ≥ χ
at the same time for all l. We focus next on such a realization.
By choosing a sufficiently large l0 , there exists η > 0 such
that
2


(46)
f0 xl+1 − f0 xl ≤ −γ l η dl , ∀l ≥ l0 .

Lemma 4. There exists a constant L̂ > 0 such that
(43)

l→∞

where in the last
inequality,
we used (45) and


= 0. Let us show
liml→∞ ∇T f0 xl − ∇T f¯0l xl
by contradiction that w.p.1. lim inf l→∞ x̄l − xl = 0.
Suppose lim inf l→∞ x̄l − xl ≥ χ > 0 with a positive

Proof: From Lemma 3, we have
0
f¯il1 (x) − f¯il2 (x) ≤ C xl1 − xl2 + e (l1 , l2 ) ,

(47)

for all x
∈
X and i
=
0, 1, ..., m, where
0
liml1 ,l2 →∞ e (l1 , l2 ) = 0. Then it follows from (47)
and (43), and the Lipschitz continuity and strong convexity
of f¯il (x) , ∀i that
0

x̄l1 − x̄l2 ≤ C1 C xl1 − xl2 + C1 e (l1 , l2 ) + C2 ν l (48)
for some constant C1 , C2 > 0. Finally, Lemma 4 follows from
(48) immediately.
Using Lemma 4 and following the same analysis as
that in [10], Proof of Theorem 1, it can be shown that
lim supl→∞ x̄l − xl = 0 w.p.1. This completes the proof.
D. Proof of Theorem 1
According to Lemma 1, Lemma 2, and (43), x∗ must be the
optimal solution of the following convex optimization problem
almost surely:

(45)

for some η > 0, where dl = x̄l − xl , liml→∞ e (l) = 0,
and the last equality follows from (44). From Assumption 1,
the gradient of f0 (x) is Lipschitz continuous, and thus there
exists L0 > 0 such that



2 l 2
f0 xl+1 ≤ f0 xl + γ l ∇T f0 xl dl + L0 γ l
d
2


2
= f0 xl + L0 γ l
dl




l
l
+ γ l ∇T f0 xl − ∇T f¯0 xl + ∇T f¯0 xl dl
2


≤ f0 xl − γ t η dl + o γ l ,

min fˆ0 (x)

(49)

x

s.t. fˆi (x) ≤ 0, i = 1, ...., m.
Since the Slater condition is satisfied, the KKT condition of
Problem (49) implies that there exist λ1 , ..., λm such that
X
∇fˆ0 (x) +
λi ∇fˆi (x∗ ) = 0,
i

fˆi (x∗ ) ≤ 0, ∀i = 1, ..., m,
λi fˆi (x∗ ) = 0, ∀i = 1, ..., m.
∗

(50)

It follows from Lemma 1 and (50) that x also satisfies the
KKT condition of Problem (7). This completes the proof.

13

R EFERENCES
[1] 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, no. 3, pp. 1499–1513, Mar. 2014.
[2] L. Liang, W. Xu, and X. Dong, “Low-complexity hybrid precoding
in massive multiuser MIMO systems,” IEEE Wireless Commun. Lett.,
vol. 3, no. 6, pp. 653–656, Dec. 2014.
[3] A. Alkhateeb, G. Leus, and R. Heath, “Limited feedback hybrid precoding for multi-user millimeter wave systems,” IEEE Trans. Wireless
Commun., vol. 14, no. 11, pp. 6481–6494, Nov. 2015.
[4] A. Liu and V. K. N. Lau, “Impact of CSI knowledge on the codebookbased hybrid beamforming in massive MIMO,” IEEE Transactions on
Signal Processing, vol. 64, no. 24, pp. 6545–6556, Dec 2016.
[5] S. Park, J. Park, A. Yazdan, and R. W. Heath, “Exploiting spatial
channel covariance for hybrid precoding in massive MIMO systems,”
IEEE Trans. Signal Processing, vol. 65, no. 14, pp. 3818–3832, July
2017.
[6] A. Liu and V. K. N. Lau, “Phase only RF precoding for massive MIMO
systems with limited RF chains,” IEEE Trans. Signal Processing, vol. 62,
no. 17, pp. 4505–4515, Sept. 2014.
[7] ——, “Two-stage constant-envelope precoding for low-cost massive
MIMO systems,” IEEE Trans. Signal Processing, vol. 64, no. 2, pp.
485–494, Jan. 2016.
[8] A. K. Sadek, W. Su, and K. J. R. Liu, “Transmit beamforming for
space-frequency coded MIMO-OFDM systems with spatial correlation
feedback,” IEEE Trans. Commun., vol. 56, no. 10, pp. 1647–1655, Oct.
2008.
[9] A. F. Molisch, V. V. Ratnam, Z. L. S. Han, S. L. H. Nguyen, L. Li,
and K. Haneda, “Hybrid beamforming for massive MIMO - a survey,”
arXiv preprint arXiv:1609.05078, 2016.
[10] Y. Yang, G. Scutari, D. P. Palomar, and M. Pesavento, “A parallel decomposition method for nonconvex stochastic multi-agent optimization
problems,” IEEE Trans. Signal Processing, vol. 64, no. 11, pp. 2949–
2964, June 2016.
[11] X. Zhang, A. Molisch, and S.-Y. Kung, “Variable-phase-shift-based RFbaseband codesign for MIMO antenna selection,” IEEE Trans. Signal
Processing, vol. 53, no. 11, pp. 4091–4103, Nov. 2005.
[12] 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,
no. 1, pp. 195 – 202, Jan. 2005.
[13] 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. Info. Theory, vol. 58, no. 7, pp.
4509–4537, Jul. 2012.
[14] A. Liu and V. Lau, “Joint power and antenna selection optimization
in large cloud radio access networks,” IEEE Trans. Signal Processing,
vol. 62, no. 5, pp. 1319–1328, March 2014.
[15] A. Adhikary, J. Nam, J.-Y. Ahn, and G. Caire, “Joint spatial division and
multiplexing - the large-scale array regime,” IEEE Trans. Info. Theory,
vol. 59, no. 10, pp. 6441–6463, Oct. 2013.
[16] F. Sohrabi and W. Yu, “Hybrid digital and analog beamforming design
for large-scale antenna arrays,” IEEE Journal of Selected Topics in
Signal Processing, vol. 10, no. 3, pp. 501–513, April 2016.
[17] D. Wipf and S. Nagarajan, “Iterative reweighted ell1 and ell2 methods
for finding sparse solutions,” IEEE Journal of Selected Topics in Signal
Processing, vol. 4, no. 2, pp. 317–329, April 2010.
[18] F. Kelly, A. Maulloo, and D. Tan, “Rate control for communication
networks: Shadow price proportional fairness and stability,” J. Oper.
Res. Soc., vol. 49, pp. 237–252, 1998.
[19] A. Liu and V. Lau, “Hierarchical interference mitigation for massive
MIMO cellular networks,” IEEE Trans. Signal Processing, vol. 62,
no. 18, pp. 4786–4797, Sept. 2014.
[20] T. V. Chien, E. BjÃ¶rnson, and E. G. Larsson, “Joint power allocation
and user association optimization for massive MIMO systems,” IEEE
Trans. Wireless Commun., vol. 15, no. 9, pp. 6384–6399, Sept 2016.
[21] S. Boyd, L. Xiao, and A. Mutapcic, “Subgradient methods,” 2003.
[Online]. Available: http://www.stanford.edu/class/ee392o
[22] S. Boyd and L. Vandenberghe, Convex Optimization.
Cambridge
University Press, 2004.
[23] M. Razaviyayn, “Successive convex approximation: Analysis and applications,” Ph.D. dissertation, University of Minnesota, 2014.
[24] A. Ruszczynski, “Feasible direction methods for stochastic programming
problems,” Math. Programm., vol. 19, no. 1, pp. 220–229, Dec. 1980.

