# Joint Precoding and RRH selection for User-centric Green MIMO C-RAN

Cunhua Pan, Huiling Zhu, Nathan J. Gomes, and Jiangzhou Wang, _Fellow_ ,
_IEEE_ This work has been accepted by IEEE TWC. Part of this work has been
presented in IEEE Globecom 2016 [1].C. Pan was with the School of Engineering
and Digital Arts, University of Kent, Canterbury, Kent, CT2 7NZ, U.K. He is
now with the School of Electronic Engineering and Computer Science, Queen Mary
University of London, London E1 4NS, U.K. (Email:{C.Pan}@qmul.ac.uk)H. Zhu, N.
Gomes and J. Wang are with the School of Engineering and Digital Arts,
University of Kent, Canterbury, Kent, CT2 7NZ, U.K. (Email:{ H.Zhu, N.J.Gomes,
J.Z.Wang}@kent.ac.uk).

###### Abstract

This paper jointly optimizes the precoding matrices and the set of active
remote radio heads (RRHs) to minimize the network power consumption (NPC) for
a user-centric cloud radio access network (C-RAN), where both the RRHs and
users have multiple antennas and each user is served by its nearby RRHs. Both
users’ rate requirements and per-RRH power constraints are considered. Due to
these conflicting constraints, this optimization problem may be infeasible. In
this paper, we propose to solve this problem in two stages. In Stage I, a low-
complexity user selection algorithm is proposed to find the largest subset of
feasible users. In Stage II, a low-complexity algorithm is proposed to solve
the optimization problem with the users selected from Stage I. Specifically,
the re-weighted $l_{1}$-norm minimization method is used to transform the
original problem with non-smooth objective function into a series of weighted
power minimization (WPM) problems, each of which can be solved by the weighted
minimum mean square error (WMMSE) method. The solution obtained by the WMMSE
method is proved to satisfy the Karush-Kuhn-Tucker (KKT) conditions of the WPM
problem. Moreover, a low-complexity algorithm based on Newton’s method and the
gradient descent method is developed to update the precoder matrices in each
iteration of the WMMSE method. Simulation results demonstrate the rapid
convergence of the proposed algorithms and the benefits of equipping multiple
antennas at the user side. Moreover, the proposed algorithm is shown to
achieve near-optimal performance in terms of NPC.

##  I Introduction

Mobile communications has been developing very rapidly [2, 3, 4]. In recent
years, C-RAN has been proposed as a promising solution to support the
exponential growth of mobile data traffic [5, 6]. In C-RAN, all the baseband
processing is performed at the baseband unit (BBU) pool with powerful
computation capacity, while the remote radio heads (RRHs) perform the basic
functionalities of signal processing [7, 8]. The RRHs are geographically
distributed away from each other, but connected to the BBU pool through
optical fiber transport links. Under the C-RAN architecture, centralized
signal processing technologies can be realized. Hence, significant performance
gains can be achieved. In addition, the RRHs can be densely deployed in the
network with low operation cost due to their simple functionalities. This will
significantly reduce the average access distance for the users, and thus
lowers the transmission power.

On the other hand, it was reported that the total energy consumption of
wireless communications contributes more than 3 percent of the worldwide
electrical energy consumption [9], and this portion is expected to grow in the
near future due to the explosive growth of high-data-rate applications and
mobile devices. Hence, energy efficiency has attracted extensive interest and
becomes one of the main performance metrics in the future fifth generation
(5G) systems [10, 11, 12, 13, 14]. When a large number of RRHs are deployed in
the network, the network power consumption (NPC) of C-RAN will become
considerable due to the increasing circuit power consumption of the RRHs.
Fortunately, it was reported in [15] that the traffic load varies
substantially over both time and space due to user mobility and varying
channel state. Hence, the NPC can be significantly reduced by putting some
RRHs with light load into sleep mode while maintaining the quality of service
(QoS) requirements of the users, which is the focus of this paper.

Recently, the NPC minimization problem for C-RAN has been extensively studied
in [16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26]. These papers formulated the
joint RRH selection and beamforming vector optimization problem as a mixed-
integer non-linear programming (MINLP) problem, which has a nonconvex
discontinuous $l_{0}$-norm in the objective function or constraints. We
summarize the existing approaches to solve the MINLP problem as follows. The
first approach was proposed in [16], which first reformulated the problem as
an extended mixed integer second-order cone programming (SOCP) and then
applied the branch-and-cut method to obtain the optimal solution. In the
second approach in [17, 18], the MINLP was first decomposed into a master
problem and a beamforming subproblem. Then, an iterative algorithm based on
the Benders decomposition was derived to find the optimal solution. Although
these two approaches yield the optimal solution, they have an exponential
complexity. The third approach is the smooth function method, where the
$l_{0}$-norm was approximated as Gaussian-like function in [19], the
exponential function in [20], and arctangent function in [21]. However, the
smooth function cannot produce sparse solutions in general. The last approach
was inspired by the compression sensing, named re-weighted $l_{1}$-norm
minimization method [27]. This method has been widely adopted in the
literature [28, 22, 23, 24, 25, 26] due to its low computational complexity
and sparsity guarantee, which will also be applied in this paper.

All of the above papers only considered the single-antenna user (SAU) case.
With the increasing development in antenna technology [29, 30], it is possible
to equip the wireless devices with multiple antennas. When both the
transmitter and the receiver are equipped with multiple antennas, multiple
streams can be transmitted simultaneously, rather than only one stream in the
SAU case. Simulation results show that with the increasing number of receive
antennas, more users can be admitted. Therefore, in this paper, we consider
the multiple-antenna user (MAU) case and jointly optimize the precoding
matrices and the set of active RRHs to minimize the NPC subject to users’ rate
requirements and per-RRH power constraints.

Unfortunately, the techniques in [16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26]
dealing with the SAU case cannot be extended directly to the MAU case. The
reasons are as follows. Firstly, since the rate constraints and power
constraints are conflicting with each other, this problem may be infeasible.
In the SAU networks, the rate requirements can be equivalently represented as
signal-to-interference-plus-noise ratio (SINR) constraints, which can be
transformed into an SOCP problem. Hence, the feasibility of the original
problem can be easily checked by solving the SOCP feasibility problem.
However, the rate constraints in the MAU case is non-convex and much more
complex due to the complicated rate expression, which cannot be transformed
into the SOCP formulation as in the SAU case. Hence, new techniques need to be
developed to check the feasibility of the original problem. Secondly, even
though the original problem is checked to be feasible, how to solve it is
still difficult, since it cannot be transformed into an SOCP problem as in the
SAU case. [31] proposed the weighted minimum mean square error (WMMSE) method
to solve the rate maximization problem for MIMO interfering broadcast
channels, where the rate expression is in the objective function. Recently,
there have been some work in applying the WMMSE method to solve the energy
efficiency (measured in bit/s/Joule) optimization problems under rate
constraints [32, 33]. However, these researches have not addressed the
feasibility problem due to the incorporated rate constraints. Only in [33], a
heuristic method was proposed to check the feasibility based on the
interference alignment technique, under the assumption that the transmit power
is approaching infinity, which in not practical. Since the problem considered
in this paper imposes power constraints at each RRH, the heuristic method
developed in [33] is not applicable. More importantly, they have not revealed
the hidden property of applying WMMSE method to the optimization problem with
rate constraints, such as the convergence property and the optimality of the
solutions.

To the best of our knowledge, this paper is the first attempt to solve the
joint RRH and precoding optimization problem to minimize the NPC in the MAU
based user-centric C-RAN, where each user can be served by an arbitrary subset
of RRHs. Due to the conflicting constraints, this problem may be infeasible.
Some users should be removed or rescheduled for the next transmission to
guarantee the rate requirements of other users. We provide a comprehensive
analysis for this problem by considering two stages: user selection in Stage I
and algorithm design in Stage II. The main contributions of this paper are
summarized as follows:

  1. 1.

In Stage I, a low-complexity user selection approach is proposed to maximize
the number of admitted users that can have their QoS requirements satisfied.
Specifically, in each step we solve an alternative problem by introducing a
series of auxiliary variables. This alternative problem is always feasible. By
replacing the rate expression in the constraints with its lower-bound, an
iterative algorithm is proposed to solve this problem along with the
complexity and convergence analysis of the algorithm. The alternative problem
should be solved at most $K$ times, where $K$ is the total number of users.
Its complexity is much lower than the optimal exhaustive user selection method
that has an exponential complexity. Simulation results show that both
algorithms achieve similar performance.

  2. 2.

In Stage II, a low-complexity algorithm is proposed to solve the NPC
minimization problem with the users selected from Stage I. Specifically, the
re-weighted $l_{1}$-norm minimization method is adopted to convert the non-
smooth optimization problem into a series of smooth weighted power
minimization (WPM) problems. We again replace the rate expression with its
lower-bound and adapt the WMMSE algorithm originally designed for a rate
maximization problem to solve the WPM problem. In addition, we strictly prove
that when the WMMSE algorithm is initialized with a feasible solution, the
sequences of precoder matrices generated in the iterative procedure will
finally converge to the Karush-Kuhn-Tucker (KKT) point of the WPM problem.

  3. 3.

In each iteration of the WMMSE algorithm, there is a subproblem for the
precoder matrices being updated with some other fixed variables. Most existing
papers [32, 33, 34, 25, 35] directly transform it into an SOCP problem and
apply the interior point method [36] to solve it, which may incur high
computational complexity. In this paper, we go one step further and develop a
low-complexity algorithm to solve this subproblem by exploiting its special
structure. Specifically, we equivalently solve its dual problem because the
subproblem is a convex problem. Fortunately, the objective function of the
dual problem is differentiable, and the block coordinate descent (BCD) method
is adopted to solve the dual problem. In each iteration of the BCD method,
Newton’s method and the gradient descent method are applied to update the
Lagrangian multipliers. It is strictly proved that the BCD method can obtain
the globally optimal solution of the subproblem. Complexity analysis in
conjunction with the simulation results show that the BCD method has a much
lower computational complexity than the interior point method.

This paper is organized as follows. In Section II, we introduce the system
model and formulate the optimization problems. In Section III, a new approach
is introduced to select the maximum number of admitted users. An iterative
algorithm with low complexity is provided in Section IV. Simulation results
are presented in Section V. Conclusions are drawn in Section VI.

Notations: Uppercase and lowercase boldface denote matrices and vectors,
respectively. For a matrix ${\bf{A}}$, ${\left\|{\bf{A}}\right\|_{F}}$ denotes
the Frobenius norm of ${\bf{A}}$ and ${{\bf{A}}^{H}}$ represents the Hermitian
transpose of ${\bf{A}}$. ${\bf{I}}_{m}$ denotes a $m\times m$ identity matrix.
For a vector ${\bf{a}}$, ${\rm{diag}}(\bf{a})$ denotes the diagonal matrix
with diagonal elements given by $\bf{a}$. ${\rm{blkdiag}}(\bf{\cdot})$
represent the block-diagonal matrices. $\mathbb{E}(\cdot)$, and
${\rm{Tr}}(\cdot)$ represent expectation, trace operators, respectively.
${\bf{A}}\succeq{\bf{B}}$ means ${\bf{A}}-{\bf{B}}$ is a positive semidefinite
matrix. For vector ${\bf{a}}\in{\mathbb{C}^{n\times 1}}$,
${\left\|{\bf{a}}\right\|_{2}}$ is the Euclidean norm. ${\cal C}{\cal
N}\left({{\bf{0}},\sigma^{2}{{\bf{I}}}}\right)$ represents the complex
circularly symmetric Gaussian distribution with zero mean vector and
covariance matrix $\sigma^{2}{{\bf{I}}}$. For a vector ${\bf{x}}$,
${\left\|{\bf{x}}\right\|_{0}}$ is ${l_{0}}$-norm, means the number of nonzero
entries in a vector.

##  II System Model and Problem Formulation

###  II-A System model

Consider a downlink C-RAN consisting of $I$ RRHs and $K$ users 111In dense
networks, the number of RRHs may be larger than the number of users so that
the average distance between serving RRHs and users can be significantly
reduced, leading to improved performance. In some extreme cases, each user may
be served by its dedicated RRHs as in [37, 38], where each RRH serves only one
user., where each RRH is equipped with $M$ transmit antennas and each user has
$N$ receive antennas, as shown in Fig. 1. Denote the set of RRHs and users as
${\cal I}=\left\\{{1,\cdots,I}\right\\}$ and $\bar{\cal
U}=\left\\{{1,\cdots,K}\right\\}$, respectively. It is assumed that each RRH
is connected to the BBU pool via fronthaul link and the BBU pool has access to
all users’ CSI and data information.

Figure 1: Illustration of a C-RAN with thirteen RRHs and six users, where
user-centric clustering technique is adopted. In this example, each user is
served by its nearby RRHs within the dotted circle centered at itself. The
RRHs that are not in any users’ candidate set are turned into idle mode, such
as RRH3 and RRH 5.

Let ${{\cal U}}\subseteq\bar{\cal U}$ be the set of users that can be admitted
to this networks. To reduce the computational complexity of the dense network,
the user-centric clustering method is adopted, where each user $k\in{\cal U}$
is assumed to be served by its nearby RRHs since the distant RRHs contribute
less to user’s signal quality due to the large path loss. The unselected RRHs
are turned into idle mode, such as RRH 3 and RRH 5 in Fig. 1. Let ${{\cal
I}_{k}}\subseteq{\cal I}$ and ${{{\cal U}}_{i}}\subseteq{{\cal U}}$ be the
candidate set of RRHs for serving user $k$ and candidate set of users served
by RRH $i$, respectively. Note that the set of RRHs serving the users may
overlap with each other. For example, in Fig. 1, RRH 12 jointly serves user 1
and user 6.

Denote ${\bf{V}}_{i,k}\in\mathbb{C}^{M\times d}$ as the precoding matrix used
by the $i$th RRH to transmit data vector ${\bf{s}}_{k}\in\mathbb{C}^{d\times
1}$ to the $k$th user, where $d$ is the number of data streams for each user,
and ${\bf{s}}_{k}$ satisfies
$\mathbb{E}\left[{{{\bf{s}}_{k}}{\bf{s}}_{k}^{\rm{H}}}\right]={\bf{I}}_{d}$
and
$\mathbb{E}\left[{{{\bf{s}}_{k}}{\bf{s}}_{l}^{\rm{H}}}\right]={\bf{0}},{\rm{for}}\
l\neq k$. Let ${{\bf{\bar{V}}}_{k}}={\left[{{\bf{V}}_{i,k}^{\rm{H}},\forall
i\in{{\cal I}_{k}}}\right]^{\rm{H}}}\in{\mathbb{C}^{\left|{{{\cal
I}_{k}}}\right|M\times d}}$ be the big precoding matrix for user $k$ from all
RRHs in ${\cal I}_{k}$. In addition, define a set of new channel matrices
${\bf{\bar{H}}}_{j,k}=[{{\bf{H}}_{i,k},\forall i\in{{\cal
I}_{j}}}]\in{{\mathbb{C}}^{N\times\left|{{{\cal I}_{l}}}\right|M}}$,
representing the overall CSI from RRHs in ${{\cal I}_{j}}$ to user $k$, where
${{{\bf{H}}_{i,k}}}\in{\mathbb{C}^{N\times M}}$ denotes the channel matrix
from the $i$th RRH to the $k$th user. Then, the received signal vector at the
$k$th user, denoted as ${\bf{y}}_{k}\in\mathbb{C}^{N\times 1}$, is given by

$\begin{array}[]{l}{\bf{y}}_{k}={{\bf{\bar{H}}}_{k,k}}{{\bf{\bar{V}}}_{k}}{{\bf{s}}_{k}}+\sum\nolimits_{j\in{\cal
U},j\neq
k}{{\bf{\bar{H}}}_{j,k}}{{\bf{\bar{V}}}_{j}}{{\bf{s}}_{j}}+{\bf{n}}_{k},\end{array}$
(1)

where ${\bf{n}}_{k}$ is the noise vector at the $k$th user, which satisfies
${\cal C}{\cal N}\left({{\bf{0}},\sigma_{k}^{2}{{\bf{I}}_{N}}}\right)$. Then,
the achievable rate (nat/s/Hz) of the $k$th user is given by[39]

$\begin{array}[]{l}{R_{k}}({\bf{V}})={\rm{lo}}{{\rm{g}}}\;\left|{{\bf{I}}+{{\bf{\bar{H}}}_{k,k}}{{\bf{\bar{V}}}_{k}}{\bf{\bar{V}}}_{k}^{H}{\bf{\bar{H}}}_{k,k}^{H}{\bf{J}}_{k}^{-1}}\right|,\end{array}$
(2)

where ${\rm{log}}(\cdot)$ is the base of natural logarithm,
${\bf{J}}_{k}=\sum\nolimits_{j\in{\cal U},j\neq
k}{{\bf{\bar{H}}}_{j,k}}{{\bf{\bar{V}}}_{j}}{{\bf{\bar{V}}}_{j}^{H}}{{\bf{\bar{H}}}_{j,k}^{H}}+{\sigma^{2}_{k}}\bf{I}$
is the interference-plus-noise covariance matrix, and ${\bf{V}}$ is the
collection of all precoding matrices. Each user’s data rate should be larger
than the minimum requirement:

${\rm{C1}:}\ {R_{k}}({\bf{V}})\geq R_{k,{\rm{min}}},\forall k\in{\cal U}.$ (3)

With densely deployed RRHs, the power consumption on the RRHs and the
corresponding fronthaul links may be significant. Switching off some RRHs and
the corresponding fronthual links may be a good option to reduce the NPC. To
this end, it is critical to model the NPC.

###  II-B NPC model

The realistic NPC model should consist of three parts: power consumption at
the RRHs, that at the fronthaul links and that at the BBU pool.

As in [22], the power consumption at RRH $i$ can be modeled as follows:

$P_{i}^{{\rm{rrh}}}({\bf{V}})=\left\\{\begin{array}[]{l}{{{\eta_{i}}}}P_{i}^{{\rm{tr}}}({\bf{V}})+MP_{i}^{{\rm{a,rrh}}},\
{\rm{if}}\ P_{i}^{{\rm{tr}}}({\bf{V}})>0\\\
MP_{i}^{{\rm{s,rrh}}},\qquad\qquad\quad\ {\rm{if}}\
P_{i}^{{\rm{tr}}}({\bf{V}})=0\end{array}\right.$ (4)

where $\eta_{i}>1$ accounts for the inefficiency of the power amplifier of RRH
$i$, $P_{i}^{{\rm{tr}}}({\bf{V}})$ is the total transmit power of RRH $i$
given by $P_{i}^{{\rm{tr}}}({\bf{V}}){\rm{=}}\sum\nolimits_{k\in{{\cal
U}_{i}}}{\left\|{{{\bf{V}}_{i,k}}}\right\|_{F}^{2}}$ that satisfies the power
constraint:

${\rm{C2}:}\ P_{i}^{{\rm{tr}}}({\bf{V}})\leq P_{i,\rm{max}},\forall i\in\cal
I,$ (5)

$P_{i}^{{\rm{a,rrh}}}$ and $P_{i}^{{\rm{s,rrh}}}$ represent the power
consumption for each antenna (or each RF chain) when RRH $i$ is in active mode
and sleep mode, respectively. In practical systems, $P_{i}^{{\rm{active}}}$ is
much higher than $P_{i}^{{\rm{sleep}}}$, which motivates us to switch off some
RRHs.

In general, more power consumption will be consumed on the fronthaul links
when they support high data rates. In [26], this power was modeled to be
proportional to the total fronthaul data rate. We modify the model in [26] to
account for the power when the fronthaul links are in the sleep mode as
follows:

$P_{i}^{{\rm{fr}}}({\bf{V}})=\left\\{\begin{array}[]{l}{\rho_{i}}\sum\nolimits_{k\in{{\cal
U}_{i}}}{R_{k}}({\bf{V}})+P_{i}^{a,{\rm{fr}}},\ {\rm{if}}\
P_{i}^{{\rm{tr}}}({\bf{V}})>0,\\\
P_{i}^{s,{\rm{fr}}},\qquad\qquad\qquad\quad\quad\ {\rm{if}}\
P_{i}^{{\rm{tr}}}({\bf{V}})=0{\rm{.}}\end{array}\right.$ (6)

where $\rho_{i}$ is the proportional factor for fronthaul link $i$. The power
consumed in the BBU pool mainly depends on the computational complexity for
signal processing. However, how to accurately model this kind of power
consumption is still not fully understood. As in most papers [16, 22, 23, 26],
the BBU power consumption is modeled as a constant ${P_{{\rm{BBU}}}}$ for
simplicity. Let ${\cal A}$ denote the active RRH set. Then, the NPC can be
modeled as

$\displaystyle\hat{P}({\cal A},{\bf{V}})$ $\displaystyle=$
$\displaystyle\sum\nolimits_{i\in{\cal
I}}{\left({P_{i}^{{\rm{rrh}}}({\bf{V}})+P_{i}^{{\rm{fr}}}({\bf{V}})}\right)}+{P_{{\rm{BBU}}}}$
(7) $\displaystyle=$ $\displaystyle\sum\nolimits_{i\in{\cal
A}}{\left({{\eta_{i}}P_{i}^{{\rm{tr}}}({\bf{V}})+{\rho_{i}}\sum\nolimits_{k\in{{\cal
U}_{i}}}{{R_{k}}({\bf{V}})}+P_{i}^{c}}\right)}+\sum\nolimits_{i\in{\cal
I}}{P_{i}^{s}}+{P_{{\rm{BBU}}}},$ (8)

where $P_{i}^{c}$ and $P_{i}^{s}$ are two constants, given by
$P_{i}^{c}=M(P_{i}^{a,{\rm{rrh}}}-P_{i}^{s,{\rm{rrh}}})+P_{i}^{a,{\rm{fr}}}-P_{i}^{s,{\rm{fr}}}$
and $P_{i}^{s}=MP_{i}^{s,{\rm{rrh}}}+P_{i}^{s,{\rm{fr}}}$.

###  II-C Problem Formulation

Due to the power constraints C2, the rate requirements C1 may not be satisfied
for all users. Some users should be removed to make the optimization problem
feasible. Hence, we formulate a two-stage optimization problem. In Stage I,
one should maximize the number of admitted users that can be supported by the
system; in Stage II, one should jointly select some RRHs and optimize the
precoding matrices to minimize the NPC with the selected users from Stage I.

Specifically, the optimization problem in Stage I can be formulated as

$\begin{array}[]{l}\mathop{\max}\limits_{{\bf{V}},{\cal
U}\subseteq{\overline{\cal U}}}\quad\left|\cal U\right|\\\
\quad{\rm{s}}.{\rm{t}}.\quad{\kern 1.0pt}{\rm{C1}},{\rm{C2}}.\end{array}$ (9)

Then in Stage II, we aim to jointly select the RRHs and optimize the precoding
matrices to minimize the NPC with the users selected from Stage I, which can
be formulated as222In general, the number of transmit antennas should be
optimized to additionally reduce the NPC as seen in the RRH power consumption
model in (4). However, the resulting problem will be much more difficult to
solve, and will be left for future work.

$\displaystyle\mathop{\min}\limits_{{\cal A},{\bf{V}}}\quad$
$\displaystyle\sum\nolimits_{i\in{\cal
A}}{\left({{\eta_{i}}P_{i}^{{\rm{tr}}}({\bf{V}})+{\rho_{i}}\sum\nolimits_{k\in{{\cal
U}_{i}^{\star}}}{{R_{k}}({\bf{V}})}+P_{i}^{c}}\right)}$ (10a) s.t.
$\displaystyle\ {\rm{C1}},\sum\nolimits_{k\in{{\cal
U}_{i}^{\star}}}{\left\|{{{\bf{V}}_{i,k}}}\right\|_{F}^{2}}\leq{P_{i,\max}},i\in{\cal
A},$ (10b) $\displaystyle\ \sum\nolimits_{k\in{{\cal
U}_{i}^{\star}}}{\left\|{{{\bf{V}}_{i,k}}}\right\|_{F}^{2}}=0,i\in{\cal
I}\backslash{\cal A},$ (10c)

where ${\cal U}_{i}^{\star}$ is the solution from Stage I. Note that when the
system parameters are given, the last two terms in (8) are constants, and are
omitted in the objective function.

Both the optimization problems in the two stages are MINLP problems and are
difficult to solve. The intuitive approach to solve this kind of problems is
through the exhaustive search. For example, to solve the NPC minimization
problem in Stage II, one must solve the precoding matrices that minimizes the
NPC with each given $\cal A$ and obtain the corresponding objective value.
Finally, the $\cal A$ that achieves the minimum NPC together with the
corresponding precoding matrices is the optimal solution of Problem (10).
However, the exhaustive search has exponentially prohibitive complexity with
respect to the number of RRHs, which is hard to be implemented in practice in
dense C-RANs. The same issue holds for the user selection problem in Stage I,
where the exhaustive search method has an exponential complexity of the number
of users. Hence, this motivates us to develop low-complexity algorithms to
solve these two Problems.

##  III Stage I: Low-complexity User Selection Algorithm

In this section, we provide a low-complexity user selection algorithm to
guarantee the rate requirements of other users. Specifically, for an arbitrary
given subset of users $\cal U$, we construct an alternative problem by
introducing a series of auxiliary variables $\\{\alpha_{k}\\}_{k\in\cal U}$:

$\begin{array}[]{*{20}{l}}{\mathop{{\rm{min}}}\limits_{{{\left\\{{{\alpha_{k}}}\right\\}}_{k\in{\cal
U}}},{\bf{V}}}\;\sum\nolimits_{k\in{\cal
U}}{{{\left({{\alpha_{k}}-1}\right)}^{2}}}}\\\
{\;{\rm{s}}.{\rm{t}}.\;\;\;{\rm{C2}},{R_{k}}({\bf{V}})\geq\alpha_{k}^{2}{R_{k,{\rm{min}}}},\forall
k\in{\cal U},}\\\ \end{array}$ (11)

Obviously, Problem (11) is always feasible and the optimal $\alpha_{k}$ for
each user $k$ should be no larger than one. This can be easily proved by
contradiction. Moreover, user $k$ can be admitted if and only if the optimal
$\alpha_{k}$ is equal to one. Hence, maximizing the number of admitted users
is equal to finding the largest subset of users $\cal U$, in which all
$\\{\alpha_{k}\\}_{k\in\cal U}$ are equal to one.

Based on the above analysis, we provide a low-complexity user selection (USC)
algorithm to solve Problem (9) in Stage I. The main idea is to remove each
user with the least $\alpha_{k}<1$ in each iteration. It is intuitive since
the user with the least $\alpha_{k}$ has the largest gap to its rate target.

Algorithm 1 USC Algorithm

1: Initialize the set of users ${\cal U}=\\{1,\cdots,K\\}$;

2: Given $\cal U$, solve Problem (11) by Algorithm 2 in Subsection III-A to
obtain $\left\\{{\alpha_{k}}\right\\}_{k\in\cal U}$ and ${{\bf{V}}}$;

3: If $\alpha_{k}=1,\forall k\in\cal U$, output ${\bf{V}}$ and ${\cal
U}^{*}\\!\\!=\\!\\!\cal U$ for the initialization of Stage II and terminate;
Otherwise, find ${k^{*}}\\!\\!=\\!\\!\arg\mathop{\min}\nolimits_{k\in\cal
U}\alpha_{k}$, remove user $k^{*}$ and update ${\cal U}={\cal
U}\backslash{k^{*}}$, go to step 2.

###  III-A Algorithm to solve Problem (11)

In step 2 of Algorithm 1, Problem (11) needs to be solved. Due to constraints
C3 in (11), Problem (11) is a non-convex problem, which is difficult to solve.
To handle this difficulty, we apply the relationships between WMMSE and the
rate expression.

We consider the linear receiver filter so that the estimated signal vector is
given by

${{\bf{\hat{s}}}_{k}}={\bf{U}}_{k}^{H}{{\bf{y}}_{k}},\forall k\in\cal U.$ (12)

where ${\bf{U}}_{k}\in\mathbb{C}^{N\times d}$ is the receiver filter of the
$k$th user. Since the signal vectors ${\bf{s}}_{k}$’s and noise
${\bf{n}}_{k}$’s are mutually independent, the mean square error (MSE) matrix
at the $k$th user is given by

$\displaystyle{{\bf{E}}_{k}}$ $\displaystyle=$
$\displaystyle{\mathbb{E}_{{\bf{s,n}}}}\left[{\left({{{{\bf{\hat{s}}}}_{k}}-{{\bf{s}}_{k}}}\right){{\left({{{{\bf{\hat{s}}}}_{k}}-{{\bf{s}}_{k}}}\right)}^{H}}}\right]$
(13) $\displaystyle=$
$\displaystyle\\!\\!\\!\\!\\!\left({{\bf{U}}_{k}^{H}{{\bf{\bar{H}}}_{k,k}}{{\bf{\bar{V}}}_{k}}\\!\\!-\\!\\!{{\bf{I}}_{d}}}\right){\left({{\bf{U}}_{k}^{H}{{\bf{\bar{H}}}_{k,k}}{{\bf{\bar{V}}}_{k}}\\!\\!-\\!\\!{{\bf{I}}_{d}}}\right)^{H}}\\!\\!+\\!\\!\sum\nolimits_{j\in{\cal
U},j\neq
k}\\!\\!{{\bf{U}}_{k}^{H}{{\bf{\bar{H}}}_{j,k}}{{\bf{\bar{V}}}_{j}}}{\bf{\bar{V}}}_{j}^{H}{\bf{\bar{H}}}_{j,k}^{H}{{\bf{U}}_{k}}\\!\\!+\\!\\!\sigma_{k}^{2}{\bf{U}}_{k}^{H}{{\bf{U}}_{k}}.$

By introducing a set of auxiliary matrices
$\\{{{\bf{W}}_{k}}\succeq{\bf{0}}\\}$, we define the following functions

${h_{k}}\left({\bf{V}},{{\bf{U}}_{k},{\bf{W}}_{k}}\right)={\log\left|{{{\bf{W}}_{k}}}\right|-{\rm{Tr}}\left({{{\bf{W}}_{k}}{{\bf{E}}_{k}}}\right)+d},\forall
k.$ (14)

where ${\bf{E}}_{k}$ is the MSE matrix of user $k$ given in (13). The
following lemma establishes the relationships between the rate expression and
function ${h_{k}}\left({\bf{V}},{{\bf{U}}_{k},{\bf{W}}_{k}}\right)$.

Lemma 1 [31] : ${h_{k}}\left({\bf{V}},{{\bf{U}}_{k},{\bf{W}}_{k}}\right)$ is a
concave function for each set of the matrices ${\bf{V}}$, ${\bf{U}}_{k}$ and
${\bf{W}}_{k}$ when the other two are given. Given ${\bf{V}}$,
${h_{k}}\left({\bf{V}},{{\bf{U}}_{k},{\bf{W}}_{k}}\right)$ is the lower-bound
of the data rate ${R_{k}}({\bf{V}})$ in (2). The optimal
${{\bf{U}}_{k},{\bf{W}}_{k}}$ for
${h_{k}}\left({\bf{V}},{{\bf{U}}_{k},{\bf{W}}_{k}}\right)$ to achieve the data
rate is given by

${\bf{U}}_{k}^{\star}={\left({\sum\nolimits_{j\in{\cal
U}}{{{\bf{\bar{H}}}_{j,k}}{\bf{\bar{V}}}_{j}{\bf{\bar{V}}}_{j}^{H}{\bf{\bar{H}}}_{j,k}^{H}}+\sigma_{k}^{2}{\bf{I}}}\right)^{-1}}{{\bf{\bar{H}}}_{k,k}}{\bf{\bar{V}}}_{k},{\bf{W}}_{k}^{\star}={\bf{E}}{{}_{k}^{\star-1}},\forall
k,\forall k$ (15)

where ${\bf{E}}_{k}^{\star}$ is obtained by plugging the expression of
${\bf{U}}_{k}^{\star}$ into the $k$th user’s MSE in (13)

${\bf{E}}_{k}^{\star}={\bf{I}}_{d}-{\bf{\bar{V}}}_{k}^{H}{\bf{\bar{H}}}_{k,k}^{H}{\left({\sum\nolimits_{j\in\cal
U}{{{\bf{\bar{H}}}_{j,k}}{\bf{\bar{V}}}_{j}{\bf{\bar{V}}}_{j}^{H}{\bf{\bar{H}}}_{j,k}^{H}}+\sigma_{k}^{2}{\bf{I}}}\right)^{-1}}{{\bf{\bar{H}}}_{k,k}}{\bf{\bar{V}}}_{k}.$
(16)

$\Box$

By replacing the first set of constraints in (11) with its lower-bound
${h_{k}}\left({\bf{V}},{{\bf{U}}_{k},{\bf{W}}_{k}}\right)$, we have the
following optimization problem

$\begin{array}[]{*{20}{l}}{\mathop{{\rm{min}}}\limits_{{{\left\\{{{\alpha_{k}}}\right\\}}_{k\in{\cal
U}}},{\bf{V}},{\bf{W}},{\bf{U}}}\;\sum\nolimits_{k\in{\cal
U}}{{{\left({{\alpha_{k}}-1}\right)}^{2}}}}\\\
{\;{\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;{\rm{C2}},{h_{k}}\left({{\bf{V}},{{\bf{U}}_{k}},{{\bf{W}}_{k}}}\right)\geq\alpha_{k}^{2}{R_{k,{\rm{min}}}},\forall
k\in{\cal U},}\end{array}$ (17)

where ${\bf{U}}$ and ${\bf{W}}$ are the collection of matrices
${{\bf{U}}_{k},\forall k}$ and ${{\bf{W}}_{k},\forall k}$, respectively.

To solve Problem (17), we apply the block coordinate descent method: given
${\bf{V}}$, update ${\bf{U}}$ and ${\bf{W}}$ by using (15); update
$\\{\alpha_{k}\\}_{k\in\cal U}$ and ${\bf{V}}$ with given ${\bf{U}}$ and
${\bf{W}}$. We only need to solve the latter one. Putting the MSE expression
in (13) into constraints C4 in Problem (17) yields

$\begin{array}[]{*{20}{l}}{\mathop{\min}\limits_{\\{\alpha_{k}\\}_{k\in\cal
U},{\bf{V}}}\;\sum\nolimits_{k\in{\cal
U}}{{{\left({{\alpha_{k}}-1}\right)}^{2}}}}\\\
\begin{array}[]{l}{\rm{s}}.{\rm{t}}.\;\
{\rm{C2}},{\rm{C5}}:{\rm{Tr}}\left({{{\left({{\bf{U}}_{k}^{H}{{\bf{\bar{H}}}_{k,k}}{{\bf{\bar{V}}}_{k}}-{{\bf{I}}_{k}}}\right)}^{H}}{{\bf{W}}_{k}}\left({{\bf{U}}_{k}^{H}{{\bf{\bar{H}}}_{k,k}}{{\bf{\bar{V}}}_{k}}-{{\bf{I}}_{k}}}\right)}\right)\\\
\qquad\quad\quad+\sum\nolimits_{j\in{\cal U},j\neq
k}{{\rm{Tr}}}\left({{\bf{\bar{V}}}_{j}^{H}{\bf{\bar{H}}}_{j,k}^{H}{{\bf{U}}_{k}}{{\bf{W}}_{k}}{\bf{U}}_{k}^{H}{{\bf{\bar{H}}}_{j,k}}{{\bf{\bar{V}}}_{j}}}\right)+{\alpha_{k}^{2}}{R_{k,\min}}\leq{t_{k}},\forall
k\in\cal U,\end{array}\end{array}$ (18)

where
${t_{k}}=\log\left|{{{\bf{W}}_{k}}}\right|+d-\sigma_{k}^{2}{\rm{Tr}}\left({{\bf{U}}_{k}^{H}{{\bf{U}}_{k}}{{\bf{W}}_{k}}}\right)$.

Without loss of generality, we assume ${\cal U}={\bar{\cal
U}}=\\{1,\cdots,K\\}$ and define the indices of ${{{\cal U}_{i}}}$ as ${{\cal
U}_{i}}=\\{q_{1}^{i},\cdots,q_{|{{\cal U}_{i}}|}^{i}\\}$. Problem (18) can be
equivalently transformed into the following problem

$\begin{array}[]{l}\mathop{\min}\limits_{\\{\alpha_{k}\\}_{k\in\cal
U},{\bf{V}}}\;\sum\nolimits_{k\in{\cal
U}}{{{\left({{\alpha_{k}}-1}\right)}^{2}}}\\\ {\rm{s}}{\rm{.t}}{\rm{.}}\
{\left\|{{{\bf{x}}_{k}}}\right\|_{2}}\leq\sqrt{{t_{k}}},\forall k\in\cal U,\\\
\quad\ \ {\left\|{{{\bf{y}}_{i}}}\right\|_{2}}\leq\sqrt{{P_{i.\max}}},\forall
i\in\cal I,\end{array}$ (19)

where ${{\bf{x}}_{k}}$ is given by

$\begin{array}[]{l}{{\bf{x}}_{k}}=\left[{{\rm{vec}}{{\left({{\bf{\bar{V}}}_{1}^{H}{\bf{\bar{H}}}_{1,k}^{H}{{\bf{U}}_{k}}{\bf{W}}_{k}^{1/2}}\right)}^{H}},\cdots,}\right.{\rm{vec}}{\left({\left({{{\bf{\bar{V}}}_{k}^{H}}{{\bf{\bar{H}}}_{k,k}^{H}}{\bf{U}}_{k}-{{\bf{I}}_{k}}}\right){\bf{W}}_{k}^{1/2}}\right)^{H}},\\\
\qquad\qquad\qquad\qquad\qquad\qquad{\left.{{\rm{}}\cdots,{\rm{vec}}{{\left({{\bf{\bar{V}}}_{K}^{H}{\bf{\bar{H}}}_{K,k}^{H}{{\bf{U}}_{k}}{\bf{W}}_{k}^{1/2}}\right)}^{H}},\alpha_{k}\sqrt{{R_{k,\min}}}}\right]^{H}}\end{array}$

and ${\bf{y}}_{i}$ is given by

${{\bf{y}}_{i}}={\left[{{\rm{vec}}{{\left({{{\bf{V}}_{i,q_{1}^{i}}}}\right)}^{H}},\cdots,{\rm{vec}}{{\left({{{\bf{V}}_{i,q_{|{{\cal
U}_{i}}|}^{i}}}}\right)}^{H}}}\right]^{H}}.$ (20)

Problem (19) is an SOCP problem for which a globally optimal solution can be
obtained by existing techniques such as interior point method [36].

Based on the above analysis, the iterative algorithm for solving Problem (11)
is formally described in Algorithm 2.

Theorem 1: Algorithm 2 will converge during the iterative procedure.

Proof: Please see Appendix A. $\Box$

Algorithm 2 Iterative Algorithm

1: Initialize iterative number $n=1$, the maximum number of iterations
$n_{\rm{max}}$. Initial precoding matrices ${\bf{V}}^{(0)}$ such that the per-
RRH power constraints are satisfied. Calculate ${\bf{U}}^{(0)}$ and
${\bf{W}}^{(0)}$ by using (15) with ${\bf{V}}^{(0)}$;

2: With ${\bf{U}}^{(n-1)}$ and ${\bf{W}}^{(n-1)}$, update
$\\{\alpha_{k}^{(n)}\\}_{k\in\cal U}$ and ${\bf{V}}^{(n)}$ by solving the SOCP
problem (19);

3: Update ${\bf{U}}^{(n)}$ and ${\bf{W}}^{(n)}$ as in (15) with
${\bf{V}}^{(n)}$;

4: If $n<n_{\rm{max}}$, set $n\leftarrow n+1$ and go to step 2. Otherwise,
terminate.

###  III-B Overall complexity to solve Problem (9) in Stage I

We first analyze the complexity of Algorithm 2 to solve Problem (11). For
simplicity, we assume that candidate size for each user is equal, i.e.,
$\left|{{{\cal I}_{k}}}\right|=l$, and ${\cal U}=\bar{\cal U}$. In each
iteration of Algorithm 2, the main complexity lies in step 2, where the SOCP
Problem (19) is solved. This problem has $2MKld+K$ real variables, $K$ SOC
constraints where each has $2Kd^{2}+1$ real dimensions, and $I$ SOC
constraints where each has $2Md\left|{{{\cal U}_{i}}}\right|$ real dimensions.
From [page 196, [40]], the complexity is
$O\left({{{\left({2MKld+K}\right)}^{2}}\left({2{K^{2}}{d^{2}}+K+2Md\sum\nolimits_{i\in{\cal
I}}{\left|{{{\cal U}_{i}}}\right|}}\right)}\right)$, and the total number of
iterations required is $O\left({\sqrt{I+K}}\right)$. Note that
$\sum\nolimits_{i\in{\cal I}}{\left|{{{\cal
U}_{i}}}\right|}{\rm{=}}\sum\nolimits_{k\in{\cal U}}{\left|{{{\cal
I}_{k}}}\right|}=Kl$, the total complexity to solve the SOCP Problem (19) is
given by
$O\left({\sqrt{I+K}{{\left({2MKld+K}\right)}^{2}}\left({2{K^{2}}{d^{2}}+K+2MdKl}\right)}\right).$
Finally, Algorithm 2 should be run at most $K$ times, then the overall
complexity to solve Problem (9) in Stage I is at most
${T_{{\rm{StageI}}}}=O\left({K\sqrt{I+K}{{\left({2MKld+K}\right)}^{2}}\left({2{K^{2}}{d^{2}}+K+2MdKl}\right)}\right).$

##  IV Stage II: A Low-complexity Algorithm to Solve Problem (10)

In this section, we provide a low-complexity algorithm to solve Problem (10)
with the selected users from Stage I. First, we adopt the re-weighted
$l_{1}$-norm method [27] to transform the original non-smooth optimization
problem into a series of WPM problems. Then, the WPM problem is solved by the
WMMSE algorithm. In each iteration of the WMMSE algorithm, there is a
subproblem that the precoder matrices should be optimized. We exploit the
special structure of the subproblem and develop a low-complexity algorithm to
solve it.

###  IV-A Reweighted $l_{1}$-norm minimization

For simplicity, the subscript in ${\cal U^{\star}}$ is omitted. It is easy to
see that the minimum rate constraints in C1 of Problem (10) hold with equality
at the optimal point, i.e., ${R_{k}}({\bf{V}})={R_{k,\min}},\forall k$. Then,
defining
$\tilde{P}_{i}^{c}\buildrel\Delta\over{=}{\rho_{i}}\sum\nolimits_{k\in{\cal
U}_{i}^{\star}}{{R_{k,\min}}}+P_{i}^{c}$ and using the ${l_{0}}$-norm, the
objective function of Problem (10) is equivalent to $\sum\nolimits_{i\in{\cal
I}}{\left({{\eta_{i}}\sum\nolimits_{k\in{{\cal
U}_{i}}}{\left\|{{{\bf{V}}_{i,k}}}\right\|_{F}^{2}+{{\left\|{\sum\nolimits_{k\in{{\cal
U}_{i}}}{\left\|{{{\bf{V}}_{i,k}}}\right\|_{F}^{2}}}\right\|}_{0}}\tilde{P}_{i}^{c}}}\right)}.$
This rewritten expression enables us to apply the compressive sensing
techniques [41], where the non-smooth ${l_{0}}$-norm objective can often be
approximated by a re-weighted ${l_{1}}$-norm, i.e.,

${\left\|{\sum\nolimits_{k\in{{\cal
U}_{i}}}{\left\|{{{\bf{V}}_{i,k}}}\right\|_{F}^{2}}}\right\|_{0}}\approx
a_{i}^{(n)}\sum\nolimits_{k\in{{\cal
U}_{i}}}{\left\|{{{\bf{V}}_{i,k}}}\right\|_{F}^{2}},\vspace{-0.2cm}$ (21)

where $a_{i}^{(n)}$ is a weight factor of the $i$th RRH at the $n$th iteration
that is iteratively updated as

$a_{i}^{(n)}=\frac{1}{{\sum\nolimits_{k\in{{\cal
U}_{i}}}{\left\|{{\bf{V}}_{i,k}^{(n)}}\right\|_{F}^{2}}+\delta}},\forall i,$
(22)

where $\delta$ is a small constant regularization parameter and
${{\bf{V}}_{i,k}^{(n)}}$ is the solution in the $n$th iteration. The above
updating rule shows that those RRHs with lower transmit power in the previous
iteration will have larger weights, which force them to be shut off in the
future iterative procedure.

By using the approximation in (21), we have the following optimization problem
that should be solved in the $n$-th iteration

$\begin{array}[]{l}\mathop{\min}\limits_{\bf{V}}\ \sum\nolimits_{i\in\cal
I}{\omega_{i}^{(n-1)}\sum\nolimits_{k\in{{\cal
U}_{i}}}{\left\|{{{\bf{V}}_{i,k}}}\right\|_{F}^{2}}}\\\
{\rm{s}}{\rm{.t.}}\quad{\rm{C1,C2}},\end{array}$ (23)

where $\omega_{i}^{(n-1)}={\eta_{i}}+{a_{i}^{(n-1)}}{\tilde{P}_{i}^{c}}.$

Based on the above analysis, the re-weighted $l_{1}$-norm based (RLN)
algorithm to solve Problem (10) is given in Algorithm 3. The convergence of
the RLN algorithm is proved in [28]. In addition, [28] showed that the RLN
algorithm is guaranteed to achieve sparse solutions, while the other smooth
approximations cannot produce sparse solutions in general.

Algorithm 3 RLN algorithm

1: Initialize a small enough $\delta$, the iterative number $n=1$, the maximum
number of iterations $N_{\rm{max}}$. Initialize ${\bf{V}}^{(0)}$ with the
outputs given by Stage I, calculate $\\{\omega_{i}^{(0)},\forall i\\}$;

2: Given $\\{\omega_{i}^{(n-1)},\forall i\\}$, solve Problem (23) to get
${\bf{V}}^{(n)}$ by using the WMMSE algorithm in Subsection IV-B;

3: Update $\\{\omega_{i}^{(n)},\forall i\\}$ with ${\bf{V}}^{(n)}$;

4: If $n\geq N_{\rm{max}}$, terminate. Otherwise, set $n\leftarrow n+1$ and go
to step 2.

###  IV-B Algorithm to Solve Problem (23)

For simplicity, the subscript of $\omega_{i}^{(n-1)}$ in Problem (23) is
omitted. The main difficulty in solving Problem (23) lies in the rate
requirement, which is non-convex. To handle this difficulty, we again apply
the relationship between WMMSE and the rate expression. Based on Lemma 1, we
replace the rate constraints in (23) with its lower bound
${h_{k}}\left({\bf{V}},{{\bf{U}}_{k},{\bf{W}}_{k}}\right)$. Define the indices
of ${{\cal I}_{k}}$ as ${{\cal
I}_{k}}=\left\\{{s_{1}^{k},\cdots,s_{\left|{\cal
I}_{k}\right|}^{k}}\right\\}$, we have the following optimization problem

$\begin{array}[]{l}\mathop{{\rm{min}}}\limits_{{\bf{V,W,U}}}\sum\nolimits_{k\in{\cal
U}}{{\rm{Tr}}\left({{\bf{\bar{V}}}_{k}^{H}}{{\bf{G}}_{k}}{{{\bf{\bar{V}}}}_{k}}\right)}\\\
\;\text{s.t.}\;\;\;{h_{k}}\left({{\bf{V}},{\bf{U}}_{k},{\bf{W}}_{k}}\right)\geq{R_{k,{\text{min}}}},\forall
k\in\cal U,\\\ \quad\quad\;\sum\nolimits_{k\in{{\cal
U}_{i}}}{\left\|{{{\bf{B}}_{i,k}}{{{\bf{\bar{V}}}}_{k}}}\right\|}_{F}^{2}\leq{P_{i,\max}},\forall
i\in{\cal I},\end{array}\vspace{-0.2cm}$ (24)

where ${{\bf{G}}_{k}}$ and ${{\bf{B}}_{i,k}}$ are both diagonal matrices,
given by

${{\bf{G}}_{k}}={\rm{blkdiag}}\left\\{{{\omega_{s_{1}^{k}}}{{\bf{I}}_{M\times
M}},\cdots,{\omega_{s_{\left|{{{\cal I}_{k}}}\right|}^{k}}}{{\bf{I}}_{M\times
M}}}\right\\}\vspace{-0.2cm}$ (25)

and

${{\bf{B}}_{i,k}}{\rm{=diag}}\left\\{{\overbrace{{{\bf{0}}_{1\times
M}}}^{s_{1}^{k}},\cdots,\overbrace{{{\bf{1}}_{1\times
M}}}^{s_{j}^{k}},\overbrace{{{\bf{0}}_{1\times
M}}}^{s_{j+1}^{k}},\cdots,\overbrace{{{\bf{0}}_{1\times M}}}^{s_{\left|{{{\cal
I}_{k}}}\right|}^{k}}}\right\\},\;{\rm{if}}\;s_{j}^{k}=i,\forall i\in{\cal
I},k\in{\cal U}.\vspace{-0.2cm}$ (26)

By solving Problem (24), we can find a solution that satisfies the KKT
conditions of Problem (23). To solve it, we again apply the block coordinate
descent method. Matrices ${\bf{U}}$ and ${\bf{W}}$ can be updated with (15).
The remaining task is to update ${\bf{V}}$ with given ${\bf{U}}$ and
${\bf{W}}$. Plugging the MSE expression in (13) into the first set of Problem
(24) yields

$\begin{array}[]{l}\mathop{{\rm{min}}}\limits_{{\bf{V}}}\
\sum\nolimits_{k\in{\cal
U}}{{\rm{Tr}}\left({{\bf{\bar{V}}}_{k}^{H}}{{\bf{G}}_{k}}{{{\bf{\bar{V}}}}_{k}}\right)}\\\
{\rm{s}}{\rm{.t}}{\rm{.}}\ \sum\nolimits_{j\in{\cal
U}}{{\rm{Tr}}\left({{\bf{\bar{V}}}_{j}^{H}{\bf{\bar{H}}}_{j,k}^{H}{{\bf{U}}_{k}}{{\bf{W}}_{k}}{\bf{U}}_{k}^{H}{{{\bf{\bar{H}}}}_{j,k}}{{{\bf{\bar{V}}}}_{j}}}\right)}\\!\\!-\\!\\!{\rm{Tr}}\left({{{\bf{W}}_{k}}{\bf{U}}_{k}^{H}{{{\bf{\bar{H}}}}_{k,k}}{{{\bf{\bar{V}}}}_{k}}}\right)\\!\\!-\\!\\!{\rm{Tr}}\left({{\bf{\bar{V}}}_{k}^{H}{\bf{\bar{H}}}_{k,k}^{H}{{\bf{U}}_{k}}{{\bf{W}}_{k}}}\right)\\!\leq\\!{c_{k}},\forall
k\\\ \quad\ \ \sum\nolimits_{k\in{{\cal
U}_{i}}}{{\rm{Tr}}}\left({{\bf{\bar{V}}}_{k}^{H}{{\bf{B}}_{i,k}}{{{\bf{\bar{V}}}}_{k}}}\right)\leq{P_{i,\max}},\forall
i,\end{array}$ (27)

where
${c_{k}}=\log\left|{{{\bf{W}}_{k}}}\right|+d-{R_{k,\min}}-{\rm{Tr}}\left({{{\bf{W}}_{k}}}\right)-\sigma_{k}^{2}{\rm{Tr}}\left({{\bf{U}}_{k}^{H}{{\bf{U}}_{k}}{{\bf{W}}_{k}}}\right)$.

Algorithm 4 WMMSE Algorithm

1: Initialize iterative number $l=1$, maximum number of iterations
$l_{\rm{max}}$, feasible ${\bf{V}}^{(0)}$, calculate ${\bf{U}}^{(0)}$ and
${\bf{W}}^{(0)}$ by using (15) with ${\bf{V}}^{(0)}$, tolerance $\varepsilon$,
calculate the objective value of Problem (24), denoted as
${\rm{Obj(}}{{\bf{V}}^{(l-1)}}{\rm{)}}$.

2: With ${\bf{U}}^{(l-1)}$ and ${\bf{W}}^{(l-1)}$, update ${\bf{V}}^{(l)}$ by
solving Problem (27) with the BCD algorithm in Subsection IV-C;

3: Update ${\bf{U}}^{(l)}$ and ${\bf{W}}^{(l)}$ as in (15) with
${\bf{V}}^{(l)}$;

4: If $l\geq l_{\rm{max}}$ or
${{\left|{{\rm{Obj(}}{{\bf{V}}^{(l-1)}}{\rm{)-Obj(}}{{\bf{V}}^{(l)}}{\rm{)}}}\right|}\mathord{\left/{\vphantom{{\left|{{\rm{Obj(}}{{\bf{V}}^{(l-1)}}{\rm{)-Obj(}}{{\bf{V}}^{(l)}}{\rm{)}}}\right|}{{\rm{Obj(}}{{\bf{V}}^{(l)}}{\rm{)}}}}}\right.\kern-1.2pt}{{\rm{Obj(}}{{\bf{V}}^{(l)}}{\rm{)}}}}<\varepsilon$,
terminate. Otherwise, set $l\leftarrow l+1$ and go to step 2.

Then, an WMMSE algorithm is proposed to solve Problem (23) in Algorithm 4. In
the following theorem, we show the property of the WMMSE algorithm.

Theorem 2: The sequence of ${\bf{V}}$ generated by the WMMSE algorithm will
converge to the KKT point of Problem (23).

Proof: Please see Appendix B. $\Box$

###  IV-C Low-complexity Algorithm to Solve Problem (27)

Since $\omega_{i}>0,\forall i$, matrices $\\{{\bf{G}}_{k},\forall k\in\cal
U\\}$ are positive definite matrices. Then, Problem (27) can be similarly
transformed an SOCP problem as in (19). Using the same method in Subsection
III-B, the total complexity to solve this problem by using the interior point
method is

${T_{{\rm{SOCP}}}}=O\left({\sqrt{I+K}{{\left({2lMKd}\right)}^{2}}\left({2{K^{2}}{d^{2}}+2dMKl}\right)}\right).$
(28)

In the following, we go one step further to design an algorithm with lower
complexity. Obviously, Problem (27) is a convex problem, and it satisfies the
Slater’s condition [36]. Hence, the duality gap between Problem (27) and its
dual problem is zero [36]. Then we can solve its dual problem instead of
directly solving it.

With some simple manipulations, the Lagrangian function of Problem (27) is
given by

$\displaystyle{\cal L}\left({{\bf{V}},{\bm{\lambda}},{\bm{\mu}}}\right)$
$\displaystyle=$ $\displaystyle\sum\nolimits_{k\in{\cal
U}}{\left({{\rm{Tr}}\left({{\bf{\bar{V}}}_{k}^{H}{{{\bf{\bar{G}}}}_{k}}{{{\bf{\bar{V}}}}_{k}}}\right)-{\rm{Tr}}\left({{\lambda_{k}}{{\bf{W}}_{k}}{\bf{U}}_{k}^{H}{{{\bf{\bar{H}}}}_{k,k}}{{{\bf{\bar{V}}}}_{k}}}\right)-{\rm{Tr}}\left({{\lambda_{k}}{\bf{\bar{V}}}_{k}^{H}{\bf{\bar{H}}}_{k,k}^{H}{{\bf{U}}_{k}}{{\bf{W}}_{k}}}\right)}\right)}$
$\displaystyle-\sum\nolimits_{k\in{\cal
U}}{{\lambda_{k}}{c_{k}}}-\sum\nolimits_{i\in{\cal
I}}{{\mu_{i}}{P_{i,\max}}},$

where ${\bm{\lambda}}=[{{\lambda_{k}},\forall k\in\cal U}]^{H}$ and
${\bm{\mu}}=[{{\mu_{i}},\forall i\in\cal I}]^{H}$ are the Lagrangian
multipliers associated with the first and second sets of constrains of Problem
(27), respectively, and ${\bf{\bar{G}}}_{k}$ is given by

${{\bf{\bar{G}}}_{k}}={{\bf{G}}}_{k}+\sum\nolimits_{j\in\cal
U}{{\lambda_{j}}{\bf{\bar{H}}}_{k,j}^{H}{{\bf{U}}_{j}}{{\bf{W}}_{j}}{\bf{U}}_{j}^{H}{{\bf{\bar{H}}}_{k,j}}}+\sum\nolimits_{i\in{\cal
I}_{k}}{{\mu_{i}}}{{\bf{B}}_{i,k}}.$

The dual function is given by

$\displaystyle g({\bm{\lambda}},{\bm{\mu}})$ $\displaystyle=$
$\displaystyle\mathop{{\rm{min}}}\limits_{\bf{V}}{\cal
L}\left({{\bf{V}},{\bm{\lambda}},{\bm{\mu}}}\right)$ (29) $\displaystyle=$
$\displaystyle\mathop{{\rm{min}}}\limits_{\bf{V}}\sum\nolimits_{k\in{\cal
U}}{\left({{\rm{Tr}}\left({{\bf{\bar{V}}}_{k}^{H}{{{\bf{\bar{G}}}}_{k}}{{{\bf{\bar{V}}}}_{k}}}\right)-{\rm{Tr}}\left({{\lambda_{k}}{{\bf{W}}_{k}}{\bf{U}}_{k}^{H}{{{\bf{\bar{H}}}}_{k,k}}{{{\bf{\bar{V}}}}_{k}}}\right)-{\rm{Tr}}\left({{\lambda_{k}}{\bf{\bar{V}}}_{k}^{H}{\bf{\bar{H}}}_{k,k}^{H}{{\bf{U}}_{k}}{{\bf{W}}_{k}}}\right)}\right)}$
$\displaystyle\qquad-\sum\nolimits_{k\in{\cal
U}}{{\lambda_{k}}{c_{k}}}-\sum\nolimits_{i\in{\cal
I}}{{\mu_{i}}{P_{i,\max}}}.$

Note that matrices $\\{{\bf{G}}_{k},\forall k\in\cal U\\}$ are positive
definite matrices. Problem (29) is a convex problem, and the optimal solution
can be obtained from its first-order derivative condition as:

${{{\bf{\bar{V}}}}_{k}}={\lambda_{k}}{\bf{\bar{G}}}_{k}^{-1}{\bf{\bar{H}}}_{k,k}^{H}{{\bf{U}}_{k}}{{\bf{W}}_{k}},\forall
k\in{\cal U}.$ (30)

By inserting this solution into (29), the dual function becomes

$g({\bm{\lambda}},{\bm{\mu}})=-\sum\nolimits_{k\in{\cal
U}}{\lambda_{k}^{2}{\rm{Tr}}\left({{\bf{W}}_{k}^{H}{\bf{U}}_{k}^{H}{{{\bf{\bar{H}}}}_{k,k}}{\bf{\bar{G}}}_{k}^{-1}{\bf{\bar{H}}}_{k,k}^{H}{{\bf{U}}_{k}}{{\bf{W}}_{k}}}\right)}-\sum\nolimits_{k\in{\cal
U}}{{\lambda_{k}}{c_{k}}}-\sum\nolimits_{i\in{\cal
I}}{{\mu_{i}}{P_{i,\max}}}.$ (31)

Hence, the dual problem of Problem (27) is given by

$\displaystyle\mathop{\max}\limits_{\\{{\lambda_{k}}\geq 0,\forall
k\\},\\{{\mu_{i}}\geq 0,\forall i\\}}g({\bm{\lambda}},{\bm{\mu}})$ (32)
$\displaystyle=$
$\displaystyle\\!\\!\mathop{\min}\limits_{\\{{\lambda_{k}}\geq 0,\forall
k\\},\\{{\mu_{i}}\geq 0,\forall i\\}}\\!\sum\nolimits_{k\in{\cal
U}}{\lambda_{k}^{2}{\rm{Tr}}\left({{\bf{W}}_{k}^{H}{\bf{U}}_{k}^{H}{{{\bf{\bar{H}}}}_{k,k}}{\bf{\bar{G}}}_{k}^{-1}{\bf{\bar{H}}}_{k,k}^{H}{{\bf{U}}_{k}}{{\bf{W}}_{k}}}\right)}\\!+\\!\sum\nolimits_{k\in{\cal
U}}{{\lambda_{k}}{c_{k}}}\\!+\\!\sum\nolimits_{i\in{\cal
I}}{{\mu_{i}}{P_{i,\max}}}$ $\displaystyle\triangleq$
$\displaystyle\mathop{\min}\limits_{\\{{\lambda_{k}}\geq 0,\forall
k\\},\\{{\mu_{i}}\geq 0,\forall i\\}}f({\bm{\lambda}},{\bm{\mu}}),$

where $f({\bm{\lambda}},{\bm{\mu}})=-g({\bm{\lambda}},{\bm{\mu}})$.

Fortunately, the objective function of the dual problem in (32) is
differentiable and the dual problem is convex [36], the descent methods such
as the gradient descent method and Newton’s method [36, 42] can be applied to
solve it. In the following, we also utilize the block coordinate descent
method to solve the dual problem (32): Optimize
$\left\\{{{\lambda_{k}}},\forall k\right\\}$ with $\left\\{{{\mu_{i}}},\forall
i\right\\}$, and vice versa.

Given $\left\\{{{\mu_{i}}},\forall i\right\\}$, Newton’s method is applied to
find the optimal $\left\\{{{\lambda_{k}}},\forall k\right\\}$ of the dual
problem, which is summarized in Algorithm 5. 333Since
$\left\\{{{\mu_{i}}},\forall i\right\\}$ are given, $f({\bm{\lambda}})$ is
short for $f({\bm{\lambda}},{\bm{\mu}})$ and the same for $f({\bm{\mu}})$
later.

Algorithm 5 Newton’s Method to Update $\left\\{{\lambda_{k}},\forall
k\right\\}$

1: Initialize iterative number $t=1$, the maximum number of iterations
$t_{\rm{max}}^{\rm{Newt}}$, initial $\bm{\lambda}^{(0)}=\bf{1}$, tolerance
$\varepsilon={10^{-10}}$, $\xi\in(0,0.5)$, $\varphi\in(0,1)$;

2: Compute the gradient $\nabla f(\bm{\lambda}^{(t-1)})$, Hessian matrix
$\nabla^{2}f(\bm{\lambda}^{(t-1)})$, the Newton direction and the decrement
$\Delta{\bm{\lambda}^{(t-1)}}=-{\left({{\nabla^{2}}f({\bm{\lambda}^{(t-1)}})}\right)^{-1}}\nabla
f({\bm{\lambda}^{(t-1)}}),o^{(t-1)}=\nabla
f{({\bm{\lambda}^{(t-1)}})^{T}}{\left({{\nabla^{2}}f({\bm{\lambda}^{(t-1)}})}\right)^{-1}}\nabla
f({\bm{\lambda}^{(t-1)}});$

3: Compute
$\bm{\bar{\lambda}}^{(t-1)}={[\bm{\lambda}^{(t-1)}+\Delta{\bm{\lambda}}^{(t-1)}]_{+}}$;

4: Update
${\bm{\lambda}}^{(t)}={\bm{\lambda}}^{(t-1)}+{\kappa^{(t-1)}}(\bm{\bar{\lambda}}^{(t-1)}-{\bm{\lambda}}^{(t-1)})$,
where ${\kappa^{(t-1)}}={\varphi^{{m^{(t-1)}}}}$ and ${m^{(t-1)}}$ is the
first non-negative integer $m$ that satisfies
$f({{\bm{\lambda}}^{(t)}})-f({{\bm{\lambda}}^{(t-1)}})\leq\xi{\varphi^{m}}\nabla
f{({{\bm{\lambda}}^{(t-1)}})^{T}}\left({{{\bar{\bm{\lambda}}}^{(t-1)}}-{{\bm{\lambda}}^{(t-1)}}}\right).$
(33)

5: If $o^{(t-1)}/2\leq\varepsilon$ or $t\geq{t_{\max}^{\rm{Newt}}}$,
terminate; Otherwise, $t\leftarrow t+1$, and go to step 2;

In step 4 of Algorithm 5, the backtracking line search method is used to find
the step size, where $\xi$ is typically chosen as a very small value and
$\varphi$ is chosen between $0$ and $1$. The step $\kappa^{(t)}$ starts with
one and then reduces by a factor of $\varphi$ until the stop condition (33) is
satisfied. Note that in each iteration of Algorithm 5, the step value
$\kappa^{(t)}$ may be different. The constant $\xi$ can be regarded as the
acceptable fraction of the decrease in the objective value of $f$ that is
predicted by the line search method.

However, to make this algorithm work, there are still problems to be solved:
how to calculate the gradient and how to compute the Hessian matrix. To derive
the expressions of the gradient and the Hessian matrix, we first introduce
some useful results in the matrix differential calculus. Given a matrix
function $\bm{\Gamma}(x)$, one has [43, 44]

$\frac{d}{{dx}}{\rm{Tr}}\left(\bm{\Gamma}(x)\right)={\rm{Tr}}\left({\frac{{d\bm{\Gamma}(x)}}{{dx}}}\right),$
(34)
$\frac{d}{{dx}}\bm{\Gamma}{(x)^{-1}}=-\bm{\Gamma}{(x)^{-1}}\frac{{d\bm{\Gamma}(x)}}{{dx}}\bm{\Gamma}{(x)^{-1}}.$
(35)

In addition, to simplify the expressions of the gradient and the Hessian
matrix, one defines some matrices:

$\begin{array}[]{l}{{{\bf{\tilde{H}}}}_{j,k}}={\bf{\bar{H}}}_{j,k}^{H}{{\bf{U}}_{k}},{{{\bf{\mathord{\buildrel{\lower
3.0pt\hbox{$\scriptscriptstyle\smile$}}\over{H}}}}}_{j,k}}={{{\bf{\tilde{H}}}}_{j,k}}{{\bf{W}}_{k}},{{{\bf{\hat{H}}}}_{j,k}}={{{\bf{\mathord{\buildrel{\lower
3.0pt\hbox{$\scriptscriptstyle\smile$}}\over{H}}}}}_{j,k}}{\bf{\tilde{H}}}_{j,k}^{H},{{{\bf{\tilde{G}}}}_{k}}={\bf{\bar{G}}}_{k}^{-1},{{\bf{C}}_{k}}={{{\bf{\tilde{G}}}}_{k}}{{{\bf{\mathord{\buildrel{\lower
3.0pt\hbox{$\scriptscriptstyle\smile$}}\over{H}}}}}_{k,k}},\\\
{{\bf{F}}_{k}}={\bf{\mathord{\buildrel{\lower
3.0pt\hbox{$\scriptscriptstyle\smile$}}\over{H}}}}_{k,k}^{H}{{\bf{C}}_{k}},{{\bf{Y}}_{j,k}}={\bf{C}}_{j}^{H}{{{\bf{\hat{H}}}}_{j,k}},{{{\bf{\tilde{Y}}}}_{j,k}}={{\bf{Y}}_{j,k}}{{{\bf{\tilde{G}}}}_{j}},{{\bf{Z}}_{j,k}}={{\bf{Y}}_{j,k}}{{\bf{C}}_{j}},\forall
j,k\in\cal U.\end{array}$ (36)

Based on the above results and definitions, the gradient can be derived as
follows:

$\nabla f(\bm{\lambda})={\left[{\frac{{\partial
f(\bm{\lambda})}}{{\partial{\lambda_{k}}}},\forall k\in\cal U}\right]^{H}},$
(37)

with

$\frac{{\partial
f(\bm{\lambda})}}{{\partial{\lambda_{k}}}}=2{\lambda_{k}}{\rm{Tr}}\left({{{\bf{F}}_{k}}}\right)-\sum\nolimits_{j\in\cal
U}{\lambda_{j}^{2}{\rm{Tr}}}\left({{{\bf{Z}}_{j,k}}}\right)+{c_{k}},k\in\cal
U.\vspace{-0.2cm}$ (38)

The Hessian matrix of $f(\bm{\lambda})$ can be calculated as:

${\left[{{\nabla^{2}}f({\bm{\lambda}})}\right]_{i,j}}=\left\\{\begin{array}[]{l}2{\rm{Tr}}({{\bf{F}}_{i}})+2\sum\nolimits_{k\in{\cal
U}}{\lambda_{k}^{2}{\rm{Tr}}\left({{{{\bf{\tilde{Y}}}}_{k,i}}{\bf{Y}}_{k,i}^{H}}\right)}-4{\lambda_{i}}{\rm{Tr}}\left({{{\bf{Z}}_{i,i}}}\right),\qquad\qquad\
\ \ {\rm{if}}\ i=j,\\\
-2{\lambda_{i}}{\rm{Tr}}\left({{{\bf{Z}}_{i,j}}}\right)-2{\lambda_{j}}{\rm{Tr}}\left({{{\bf{Z}}_{j,i}}}\right)+2\sum\nolimits_{k\in{\cal
U}}{\lambda_{k}^{2}{\mathop{\rm
Re}\nolimits}\left\\{{{\rm{Tr}}\left({{{{\bf{\tilde{Y}}}}_{k,j}}{\bf{Y}}_{k,i}^{H}}\right)}\right\\}}{\rm{,if}}\;j>i,\\\
{\left[{{\nabla^{2}}f({\bm{\lambda}})}\right]_{j,i}},\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\
\qquad\ \quad\ {\rm{if}}\ j<i.\end{array}\right.$ (39)

Next, given $\left\\{{{\lambda_{k}}},\forall k\in\cal U\right\\}$, we solve
the dual problem (32) to update $\left\\{{{\mu_{i}}},\forall i\in\cal
I\right\\}$. Here, the gradient descent method [36] is applied. Although
Newton’s method converges faster than the gradient descent method, simulation
results show that the gradient method also converges within five iterations
but it has much lower computational complexity than Newton’s method since it
does not require the calculations of the Hessian matrix and the inverse of the
Hessian matrix. The gradient descent method to update
$\left\\{{{\mu_{i}}},\forall i\in\cal I\right\\}$ is given in Algorithm 6.

In Algorithm 6, the gradient $\nabla f(\bm{\mu})$ is required. Define
${{\bf{D}}_{k}}={{\bf{C}}_{k}}{\bf{C}}_{k}^{\rm{H}}$. Then, by using the
results in (34) and (35), the gradient $\nabla f(\bm{\mu})$ can be calculated
as

$\nabla f(\bm{\mu})={\left[{\frac{{df(\bm{\mu})}}{{d{\mu_{i}}}},\forall
i\in\cal I}\right]^{H}},\vspace{-0.2cm}$ (40)

with

$\frac{{df(\bm{\mu})}}{{d{\mu_{i}}}}=-\sum\nolimits_{k\in{{\cal
U}_{i}}}{\lambda_{k}^{2}{\rm{Tr}}\left({{{\bf{B}}_{i,k}}{{\bf{D}}_{k}}}\right)}+{P_{i,\max}},\forall
i\in{\cal I}.$ (41)

Algorithm 6 Gradient Descent Method to Update
$\left\\{{{\mu_{i}}}\right\\}_{i=1}^{I}$

1: Initialize iterative number $t=1$, maximum number of iterations
$t_{\rm{max}}^{\rm{Grad}}$, initial $\bm{\mu}^{(0)}=\bf{1}$, accuracy
$\varepsilon$;

2: Compute the gradient $\nabla f(\bm{\mu}^{(t-1)})$;

3: Compute $\bm{\bar{\mu}}^{(t-1)}={[\bm{\mu}^{(t-1)}-\nabla
f(\bm{\mu}^{(t-1)})]_{+}}$;

4: Update
${\bm{\mu}}^{(t)}={\bm{\mu}}^{(t-1)}+{\kappa^{(t-1)}}(\bm{\bar{\mu}}-{\bm{\mu}}^{(t-1)})$,
where ${\kappa^{(t-1)}}={\beta^{{l^{(t-1)}}}}$ and ${l^{(t-1)}}$ is the first
non-negative integer $l$ that satisfies
$f({{\bm{\mu}}^{(t)}})-f({{\bm{\mu}}^{(t-1)}})\leq\delta{\beta^{l}}\nabla
f{({{\bm{\mu}}^{(t-1)}})^{T}}\left({{{\bar{\bm{\mu}}}^{(t-1)}}-{{\bm{\mu}}^{(t-1)}}}\right).$

5: If $t\geq{t_{\max}}$ or
${{\left|{f({\bm{\mu}^{(t)}})-f({\bm{\mu}^{(t-1)}})}\right|}\mathord{\left/{\vphantom{{\left|{f({\mu^{(t)}})-f({\mu^{(t-1)}})}\right|}{\left|{f({\mu^{(t)}})}\right|}}}\right.\kern-1.2pt}{\left|{f({\bm{\mu}^{(t)}})}\right|}}<\varepsilon$,
stop; Otherwise, $t\leftarrow t+1$, and go to step 2;

Finally, based on the above analysis, the method to solve the dual problem
(32) is given in Algorithm 7, which is named as Block Coordinate Descent (BCD)
method.

Algorithm 7 BCD Method to Solve the Dual Problem (32)

1: Initialize iterative number $n=1$, the maximum number of iterations
$n_{\rm{max}}$, initial $\bm{\lambda}^{(0)}=\bf{1}$ and
$\bm{\mu}^{(0)}=\bf{1}$, error tolerance $\varepsilon$;

2: Given $\bm{\mu}^{(n-1)}$, apply Newton’s method in Algorithm 5 to update
$\bm{\lambda}^{(n)}$;

3: Given $\bm{\lambda}^{(n)}$, employ the gradient descent method in Algorithm
6 to update $\bm{\mu}^{(n)}$;

4: If $n\geq{n_{\max}}$ or
${{\left|{f({\bm{\lambda}^{(n)}},{\bm{\mu}^{(n)}})-f({\bm{\lambda}^{(n-1)}},{\bm{\mu}^{(n-1)}})}\right|}\mathord{\left/{\vphantom{{\left|{f({\lambda^{(n)}},{\mu^{(n)}})-f({\lambda^{(t)}},{\mu^{(n-1)}})}\right|}{\left|{f({\lambda^{(n)}},{\mu^{(n)}})}\right|}}}\right.\kern-1.2pt}{\left|{f({\bm{\lambda}^{(n)}},{\bm{\mu}^{(n)}})}\right|}}<\varepsilon$,
terminate; Otherwise, set $n\leftarrow n+1$, and go to step 2;

Theorem 3: The sequences of $\bm{\mu}$ and $\bm{\lambda}$ generated by
Algorithm 7 will converge to the globally optimal solution of the dual problem
(32).

Proof: Please see Appendix C. $\Box$

When the optimal $\bm{\mu}$ and $\bm{\lambda}$ are obtained by using Algorithm
7, the optimal solution to Problem (27) is given by (30). As there is zero
duality gap between the primal problem (27) and dual problem (32), which means
that this solution is the globally optimal solution of Problem (27).

###  IV-D Overall Complexity to Solve Problem (27) in Stage II

In this subsection, we analyze the overall complexity to solve Problem (27).
It mainly includes three layers of iterations: the first layer is the RLN
algorithm to deal with the non-smooth $l_{0}$ norm, the second layer is the
WMMSE algorithm to deal with the non-convex rate constraints, and the third
layer is the BCD algorithm to solve Problem (27).

We first analyze the complexity of the third layer for BCD algorithm. Note
that ${{{\bf{\tilde{H}}}}_{j,k}}$, ${{{\bf{\mathord{\buildrel\lower
3.0pt\hbox{$\scriptscriptstyle\smile$}\over{H}}}}}_{j,k}}$, and
${{{\bf{\hat{H}}}}_{j,k}}$ can be calculated before the iterations of the BCD
Algorithm. The main complexity of the BCD Algorithm lies in step 2 and step 3,
where Newton’s method and gradient descent method are used to update
${\bm{\lambda}}$ and $\bm{\mu}$, respectively.

We first analyze the computational complexity of Newton’s method under the
same assumption in Subsection III-B. The main complexity in each iteration of
Newton’s method lies in step 2 and step 4 of Algorithm 5. We first analyze
step 2 of Algorithm 5. According to [45], the complexity of calculating
$\\{{{\bf{\tilde{G}}}}_{k},\forall k\in\cal U\\}$ is on the order of
$O\left({{{K\left({Ml}\right)}^{2.376}}}\right)$. For any two matrices
${\bf{X}}\in{{\mathbb{C}}^{m\times n}},{\bf{Y}}\in{{\mathbb{C}}^{n\times p}}$,
the complexity of computing ${\bf{XY}}$ is on the order of
$O\left({mnp}\right)$ [36]. In general, $d\ll MI$. Then, the total complexity
of computing $\left\\{{{{\bf{C}}_{k}},{{\bf{F}}_{k}},\forall k\in\cal
U}\right\\}$ is on the order of $O\left({K{M^{2}}{l^{2}}d}\right)$. Similarly,
the total complexity of computing
$\left\\{{{{\bf{Y}}_{j,k}},{{{\bf{\tilde{Y}}}}_{j,k}},{{\bf{Z}}_{j,k}}}\forall
j,k\in\cal U\right\\}$ is on the order of
$O\left({{K^{2}}{M^{2}}{l^{2}}d}\right)$. Hence, the total complexity of
computing $\left\\{{{{\bf{C}}_{k}},{{\bf{F}}_{k}}},\forall k\right\\}$ and
$\left\\{{{{\bf{Y}}_{j,k}},{{{\bf{\tilde{Y}}}}_{j,k}},{{\bf{Z}}_{j,k}},\forall
j,k\in\cal U}\right\\}$ is on the order of
$O\left({{K^{2}}{M^{2}}{l^{2}}d}\right)$. With a similar analysis, the total
complexity of computing (39) is on the order of
$O\left({{K^{3}}Ml{d^{2}}}\right)$. In addition, the complexity of computing
the inverse of $\nabla^{2}f(\bm{\lambda}^{(t)})$ is on the order of
$O\left({{K^{2.376}}}\right)$ [45]. Hence, the total complexity of step 2 of
Algorithm 5 is on the order of
$O\left({\max\left\\{{{K^{3}}Ml{d^{2}},{{\left(K{Ml}\right)}^{2.376}},{K^{2}}{{\left({Ml}\right)}^{2}}d}\right\\}}\right)$.
In the $t$th iteration of step 4 of Algorithm 5, $f({{\bm{\lambda}}^{(t+1)}})$
is required to calculate $m^{(t)}$ times. The complexity in each time is on
the order of
$O\left({\max\left\\{{{{\left({Ml}\right)}^{2.376}},K{{\left({Ml}\right)}^{2}}d}\right\\}}\right)$.
Thus, the total complexity of step 4 of Algorithm 5 is on the order of
$O\left({m^{(t)}}{\max\left\\{{{{\left({Ml}\right)}^{2.376}},K{{\left({Ml}\right)}^{2}}d}\right\\}}\right)$.
Simulation results show that in general $m^{(t)}$ is always equal to one,
which means that $f({{\bm{\lambda}}^{(t+1)}})$ only needs to be computed for
once. Hence, the complexity of step 4 of Algorithm 5 can be approximately by
$O\left({\max\left\\{{{{\left({Ml}\right)}^{2.376}},K{{\left({Ml}\right)}^{2}}d}\right\\}}\right)$.
As a result, the total complexity of Newton’s method is

$T_{\rm{Newton}}=O\left(t_{\rm{max}}^{\rm{Newt}}{\max\left\\{{{K^{3}}Ml{d^{2}},{{\left({Ml}\right)}^{2.376}},{K^{2}}{{\left({Ml}\right)}^{2}}d}\right\\}}\right).\vspace{-0.2cm}$
(42)

Simulation results show that Newton’s method converges very rapidly and in
general five iterations are enough for the algorithm to converge.

By using the similarly analytical technique to Newton’s method, the total
complexity of the gradient descent method is given by

$T_{\rm{Grad}}=O\left({t_{{\rm{max}}}^{{\rm{Grad}}}\max\left\\{{{{\left({MI}\right)}^{2.376}},K{{\left({MI}\right)}^{2}}d}\right\\}}\right).$
(43)

The simulation results in the next section show that the gradient descent
method usually converges within five iterations. Hence, in each iteration of
the BCD Algorithm, the complexity of Newton’s method dominates the complexity
of the gradient descent method.

Based on the above analysis, the overall complexity to solve Problem (27) in
Stage II is

${T_{{\rm{StageII}}}}={t_{{\rm{RLN}}}}{t_{{\rm{WMMSE}}}}{t_{{\rm{BCD}}}}\left({{T_{{\rm{Newton}}}}+{T_{{\rm{Grad}}}}}\right),$
(44)

where $t_{{\rm{RLN}}}$, $t_{{\rm{WMMSE}}}$ and $t_{{\rm{BCD}}}$ represent the
average number of iterations required by the RLN, WMMSE, and BCD algorithms,
respectively. Simulation results show that these three algorithms converge
very fast and generally five iterations are enough to achieve large portion of
the final performance.

##  V Simulation Results

In this section, we present simulation results to evaluate the performance of
the proposed algorithms. To be more realistic, we consider a wrap-around
system model shown in Fig. 3 as in [46], where the C-RAN network is deployed
in the central square with $[-1000\ {\rm{1000}}]\times[-1000\ {\rm{1000}}]$
meters, surrounded by eight uncoordinated square macrocells. It is assumed
that all the users and RRHs are uniformly and independently distributed in the
C-RAN region. We adopt the channel model that consists of four parts: 1) the
long term evolution (LTE) standard path loss model:
$PL_{i,k}=148.1+37.6{\log_{10}}d_{i,k}\ ({\rm{dB}})$, where $d_{i,k}$ (in km)
is the distance from the $i$th RRH to the $k$th user; 2) Log-normal shadowing
with zero mean and 8 dB standard derivation; 3) Rayleigh fading with zero mean
and unit variance; 4) transmit antenna power gain of 9 dBi. Each user is
assumed to have the same rate requirement, i.e.,
$R_{\rm{min}}=R_{k,\rm{min}},\forall k$, and each RRH has the same power
constraint, i.e., $P_{\rm{max}}=P_{i,\rm{max}}=4{\rm{W}},\forall i\in\cal I$.
It is assumed that each user is potentially served by its nearest $X$ RRHs,
i.e., $|{{\cal I}_{k}}|=X,\forall i$. Unless stated otherwise, the system
parameters are set as follows: error tolerance is $\varepsilon=10^{-3}$,
thermal noise power is $\sigma^{2}=-104\ {\rm{dBm}}$, $I=12$, $K=8$, $X=3$,
$M=2$, $N=2$, $d={\rm{min}}\\{M,N\\}$, $\eta_{i}=4$ [47], $\rho_{i}=0.5$ [26],
$P_{i}^{a,{\rm{rrh}}}=3.4{\rm{W}}$, $P_{i}^{s,{\rm{rrh}}}=2.15{\rm{W}}$,
$P_{i}^{a,{\rm{fr}}}=3.85{\rm{W}}$, $P_{i}^{s,{\rm{fr}}}=0.75{\rm{W}}$,
${P_{{\rm{BBU}}}}=20W$ [48, 22]. Moreover, let $\cal L$ be the set of
uncoordinated base stations (BSs) in C-RAN’s nearby eight macrocells. The
noise power at user $k$ can be modeled as
$\sigma_{k}^{2}={\sigma^{2}}+\sum\nolimits_{m\in{\cal
L}}{{P_{\max}}P{L_{m,k}}{S_{m,k}}{G_{m}}}$ [46], where $P{L_{m,k}}$ and
$S_{m,k}$ are the large-scale fading and shadowing respectively from the BS in
macrocell $m$ to user $k$, ${G_{m}}$ represents the antenna gain.

Figure 2: Illustration of a wrap-round C-RAN system model, where C-RAN is
deployed in the center of the region, which is surrounded by eight nearby
cells.

Figure 3: Convergence behaviour of Algorithm 2.

###  V-A Properties of the Proposed Algorithms

####  V-A1 Convergence behavior of Algorithm 2

Fig. 3 shows the convergence behaviour of Algorithm 2 for different numbers of
receive antennas. The results are obtained by averaging over 100 channel
realizations. Due to the non-convexity of Problem (11), different initial
points for Algorithm 2 may yield different solutions. To investigate this
effect, we consider two initialization schemes: 1) SVD-initial, in which the
beam directions for each user are chosen as the unitary matrices obtained by
the singular value decomposition (SVD) of channel matrices and the total power
at each RRH is equally allocated to the users potentially served by each RRH;
2) Rand-initial, in which both the beam directions and power allocations are
randomly generated. It can be seen from Fig. 3 that the objective value of
Problem (11) monotonically decreases during the iterative procedure for two
initialization schemes. In addition, the algorithm converges very fast and in
general six iterations are sufficient to achieve a large proportion of the
converged value for different numbers of receive antennas and different
initialization schemes. It is interesting to find that the algorithm under two
different initialization schemes will converge to almost the same value. As
expected, the converged objective value decreases with the number of receive
antennas since more degrees of freedom are available.

Figure 4: Average number of admitted users versus the rate requirements.

Figure 5: (a) Total power consumption versus the number of iterations; (b) The
number of active RRHs versus the number of iterations, where $R_{\rm{min}}=2$
nats/s/Hz.

####  V-A2 User selection performance of USC algorithm

Fig. 5 compares the performance of the USC algorithm with two algorithms:
greedy search method and exhaustive search method. For the greedy search
method, in each time we compute the objective value of Problem (11) when
excluding one user, then the user yielding the smallest objective value will
be removed. This procedure continues until all remaining users are feasible.
Note that this algorithm increases quadratically with $K$. The exhaustive
search method checks all feasible sets of users and chooses the largest one.
Its complexity increases exponentially with $K$. As expected, the number of
admitted users decreases with the rate requirements for all algorithms. The
greedy search method achieves almost the same performance as the exhaustive
one, and the performance gap between the exhaustive search algorithm and the
proposed USC algorithm can be negligible. However, the complexity of our
proposed USC algorithm only increases linearly with $K$. The impact of initial
points is also studied and we find that both initialization schemes (SVD-
initial and rand-initial) have similar performance, which is not shown here
for clarity.

####  V-A3 Convergence behaviour of the RLN algorithm

Figure 6: The converged state of one randomly generated system configuration.
The boundary user 8 is not selected as it is far from the RRHs.

Figure 7: The converged state of one randomly generated system configuration.
The boundary user 8 is not selected as it is far from the RRHs.

The convergence behaviours of the RLN algorithm are shown in Figs. 5 (a) and
(b) for the NPC and the number of the remaining RRHs in each iteration,
respectively. Three different values of $\delta$ are tested, i.e.,
$\delta={10^{-4}},{10^{-5}}$ and ${10^{-6}}$. One randomly generated channel
is used to obtain the convergence behaviour, where the USC algorithm is first
executed to find the largest feasible set of users. In this example, User 8 is
removed to guarantee the feasibility of the other users as seen in Fig. 7. It
can be seen from the figures that for all values of $\delta$, both the number
of active RRHs and the NPC decrease rapidly and there is no additional
decrease after the fifth iteration. At the converged state, only six RRHs are
active. Compared to the full cooperation strategy where all RRHs are active,
we can save large amount of power as seen from Fig. 5 (a). Fig. 7 illustrates
the converged state of the system. It can be seen that RRH 2 is switched off
since it is far from the users and User 8 is not selected as it is far from
the RRHs. We also study the impact of initialization schemes on the
performance of the RLN algorithm. The initial precoders for the RLN algorithm
are the outputs of the USC algorithm which is initialized with the SVD-initial
and rand-initial schemes. The simulation results show they achieve almost the
same performance, which is not shown here for clarity.

####  V-A4 Convergence behaviour of the WMMSE algorithm

Figure 8: Convergence behaviour of the BCD algorithm for the first iteration
of the WMMSE algorithm.

Figure 9: (a) Convergence behaviour of Newton’s method; (b) Convergence
behaviour of gradient descent method.

In step 2 of each iteration of the RLN algorithm, we need to solve Problem
(23) by using the WMMSE algorithm. Fig. 7 shows the convergence performance of
the WMMSE algorithm for the first three iterations of the RLN algorithm. It is
observed that the WMMSE algorithm converges within ten iterations for the
first iteration of RLN algorithm. However, the objective values stay almost
fixed for the second and third iterations of the RLN algorithm. This means
that only in the first iteration of RLN algorithm, some iterations are
required for the WMMSE algorithm.

####  V-A5 Convergence behaviour of the BCD algorithm

In step 2 of each iteration of the WMMSE algorithm, Problem (27) should be
solved to update the precoding matrices by using the BCD algorithm. Fig. 9
shows the convergence behaviour of the BCD algorithm for the first iteration
of the WMMSE algorithm. It is seen that the algorithm converges very fast and
one iteration is sufficient to achieve a large portion of the converged value
(99.2% in this example).

####  V-A6 Convergence behaviour of Newton’s method and the gradient descent
method

In each iteration of the BCD algorithm, Newton’s method is required to update
$\left\\{{\lambda_{k}},\forall k\right\\}$ and the gradient descent method is
applied to update $\left\\{{\mu_{i}},\forall i\right\\}$. The convergence
behaviours of these two algorithms for the first three iterations of the BCD
algorithm are shown in Figs. 9 (a) and (b), respectively. Newton’s method
requires several iterations to converge only in the first iteration of the BCD
algorithm, while stays almost constant for the second and third iterations of
the BCD algorithm. Interestingly, the gradient descent method only requires
one iteration to converge in the first iteration of the BCD algorithm and
keeps fixed during the rest of the iterations of the BCD algorithm. By
combining the complexity analysis in (42), (43) and the above convergence
behaviours, we can conclude that the BCD algorithm has a much lower
computational complexity than directly solving the SOCP problem.

####  V-A7 Impacts of the number of data streams

Figure 10: The number of admitted users versus rate requirements for different
numbers of data streams with $M=N=4$.

Figure 11: (a) NPC versus the rate requirements; (b) The corresponding number
of active RRHs versus the rate requirements.

In Fig. 11, the impact of the number of data streams on the number of admitted
users is studied. As expected, the number of admitted users decreases with the
rate requirements and larger number of data streams can support more users. We
find significant performance gains can be achieved when the number of data
streams increases from 1 to 2, especially for the high data rate requirements.
However, only marginal performance gains are achieved by the case of $d=4$
over the case of $d=2$, which comes at the higher cost of computational
complexity. This reveals that the performance saturates with the increase of
data streams $d$. In Fig. 11, the impacts of data streams on the NPC and on
the number of active RRHs are studied with the same setup in Fig. 11. For fair
comparison, we only consider the set of users that can be supported under the
case of $d=1$, so that all cases can support the selected users. Fig. 11 (a)
shows that the NPC first increases with the rate requirements when
$R_{\rm{min}}\leq 3\ {\rm{nats/s/Hz}}$ and then decreases significantly when
$R_{\rm{min}}>3\ {\rm{nats/s/Hz}}$. The reason can be explained as follows.
When $R_{\rm{min}}$ increases from $1$ to $3\ {\rm{nats/s/Hz}}$, the number of
admitted users almost keeps stable as shown in Fig. 11, while the fronthaul
power increases when the rate requirement increases and the number of active
RRHs increases to support the higher rate requirements as seen in Fig. 11 (b),
which in turn consumes more power consumption. On the other hand, when
$R_{\rm{min}}$ increase from $3$ to $6\ {\rm{nats/s/Hz}}$, the number of
admitted users decreases dramatically as shown in Fig. 11, which leads to
reduced transmit power and a reduced number of active RRHs as shown in Fig. 11
(b). Again, it is observed from Fig. 11 (a) that a greater number of data
streams requires lower NPC, but the performance gain shrinks with the number
of data streams.

####  V-A8 Impacts of the number of transmit antennas

Figure 12: The number of admitted users versus rate requirements for different
numbers of transmit antennas with $N=2$.

Figure 13: (a) NPC versus the rate requirements; (b) The corresponding number
of active RRHs versus the rate requirements.

In Fig. 13, the impact of the number of transmit antennas on the number of
admitted users is studied. As expected, the number of admitted users increases
with the number of transmit antennas due to more degrees of freedom.
Significant performance gains can be achieved by the case of $M=2$ over the
case of $M=1$, especially in the high rate regime. However, the performance
gain shrinks for the case of $M=4$ over the case of $M=2$. In Fig. 13, the
impacts of the number of antennas on the NPC and the number of active RRHs are
investigated. For fair comparison, it is also assumed that the set of users
selected from the case of $M=1$ are the input of Stage II for all cases of
different values of $M$ so that the selected users are the same and feasible
for all cases. It is interesting to find that when $M$ increases, the NPC
increases while the number of active RRHs decreases. This is mainly due to the
fact that the RRH power consumption model in (4) increases linearly with $M$,
and this increased power consumption dominates the reduced power consumption
resulting from the reduced number of active RRHs. It should be emphasized that
in some other cases with different values of system parameters, the NPC may
not increase with $M$ and the counter part happens, such as the case of the
low circuit power consumption for each antenna and high power consumption
associated with the fronthaul power consumption.

####  V-A9 Impacts of the candidate size

Figure 14: The number of admitted users versus rate requirements for different
candidate sizes.

Figure 15: (a) NPC versus the rate requirements; (b) The corresponding number
of active RRHs versus the rate requirements.

The impact of candidate size on the number of admitted users is illustrated in
Fig. 15 for a dense network with 20 RRHs and 12 users. As expected, larger
candidate sizes can support more users due to the increased degrees of
freedom. However, the performance gains decreases with the candidate sizes,
which implies that there is no need to consider distant RRHs for each user
since they contribute less to their signal strength. In general, the candidate
size should be no larger than 4 to achieve a good tradeoff between performance
and complexity. Similarly to the trend observed in Fig. 11 (a), it is seen
from Fig. 15 (a) that the NPC increases in the low rate regime, while
decreasing significantly in the high rate regime. For the former part, the
reason is that the increased fronthaul power dominates the reduced circuit
power for the reduced active RRHs. While for the latter part, the reason is
the opposite. Also, it is observed that the NPC performance gain for larger
candidate size is more obvious in the low rate regime, while the performance
is almost the same in the high rate regime. This is mainly due to the fact
that in the high rate regime, only a small number of users can be admitted,
and these users are separated far away. As a result, the multiuser
interference is not so significant and each user’s nearest RRH is able to
serve it with the rate requirement.

###  V-B Performance comparison

We compare the performance of the RLN algorithm with the following RRH
selection methods:

  * •

Exhaustive search (Exhau-search) method: For each given active RRH set ${\cal
A}$, this method first checks its feasibility. If feasible, the method will
use the WMMSE algorithm to solve the corresponding transmit power minimization
problem. The complexity of this method increase exponentially with $I$, which
is served as the performance benchmark for our proposed algorithm.

  * •

Successive RRH selection (Succesive-sel) method: This method first lets all
the RRHs be active and check its feasibility. If feasible, the method applies
the WMMSE algorithm to solve the transmit power minimization problem. Then,
the method gradually removes the RRHs according to their transmit power from
the lowest to the highest until the problem becomes infeasible. The complexity
of this scheme increases linearly with $I$.

  * •

Greedy search method: In each step, we exclude each RRH and calculate the NPC
when the remaining RRHs are active. Then, we remove the RRH so that the
remaining RRHs yield the least NPC. This procedure terminates until the
problem becomes infeasible. The complexity of this scheme increases
quadratically with $I$.

  * •

Full cooperative (Full-coop) method: In this method, all the selected RRHs in
cluster-formation stage are active and the WMMSE algorithm is used to solve
the transmit power minimization problem.

For fair comparison, we assume in the following simulation results, only the
channel realizations that are feasible for all users are considered.

####  V-B1 Impact of the rate requirements

Figure 16: (a) NPC versus the rate requirements; (b) The corresponding average
number of active RRHs versus the rate requirements. The candidate size is
$X=4$.

Figure 17: (a) NPC versus the number of RRHs; (b) The corresponding number of
active RRHs versus the number of RRHs with $R_{\rm{min}}=3\ {\rm{nats/s/Hz}}$
and $X=4$.

Figs. 17 (a) and (b) illustrate the average NPC and the corresponding number
of active RRHs versus the rate requirements, respectively. Fig. 17 (a) shows
that the RLN algorithm outperforms the ‘Succesive-sel’ method and ‘Full-coop’
method for all rate regimes. However, the performance of the ‘Greedy search’
method is slightly better than the RLN algorithm when $R_{\rm{min}}\leq
3{\rm{nats/s/Hz}}$, while the RLN algorithm outperforms the ‘Greedy search’
method in the high rate regime and the performance gain increases with the
rate requirements. Fig. 17 (b) shows a similar trend in terms of the number of
active RRHs. Compared with the optimal ‘Exhau-search’ method, the performance
loss in power consumption is at most $8\%$ when $R_{\rm{min}}=1$ nats/s/Hz,
and this gap gradually diminishes with the increase of rate requirements. In
particular, the performance gain provided by the ‘Exhau-search’ method over
the RLN algorithm is negligible when $R_{\rm{min}}=5$ nats/s/Hz. As expected,
the ‘Full-coop’ method consumes the highest power since all selected RRHs are
active.

####  V-B2 Impact of the number of RRHs

Figs. 17 (a) and (b) illustrate the average NPC and the corresponding number
of active RRHs versus the total number of RRHs, respectively. It is seen that
the NPC achieved by all schemes decreases with $I$ due to the fact that when
there are more RRHs, the average access distance between users and RRHs
decreases significantly and thus leads to more reduced transmit power. It is
again observed that the performance of the RLN algorithm is superior to that
of the ‘Succesive-sel’ method. This implies that selecting the RRHs only based
on the transmit power is not enough, and may incur significant performance
loss. Note that the ‘Greedy search’ method requires higher power consumption
than the RLN algorithm for all numbers of RRHs, especially when $I=6$. Also,
the performance of ‘Exhau-search’ method is slightly better than the RLN
algorithm. Note that although the number of active RRHs increases slightly
with the total number of RRHs as seen in Figs. 17 (b), the NPC decreases. This
may due to the fact that the overall transmit power reduction overwhelms the
increase of circuit power.

##  VI Conclusion

In this paper, a joint selection of active RRHs and optimization of the
precoding matrices which minimizes the NPC for the MIMO C-RAN, while
guaranteeing users’ rate requirements and per-RRH power constraints, has been
studied. A low-complexity user selection was proposed to guarantee the
feasibility of the other users. Then a low-complexity iterative algorithm,
based on the reweighted $l_{1}$-norm minimization method, WMMSE algorithm,
Newton’s method, and gradient descent method, was proposed to solve the
network power minimization problem for the selected users. Simulation results
show that the proposed algorithms converge fast, which is attractive for
practical implementation. Also, more antennas at the user side can admit more
users. The proposed user selection algorithm was shown to achieve the similar
performance as the optimal exhaustive search method. Moreover, our proposed
algorithm was shown to achieve much greater power savings than the full
cooperation method, and the performance loss compared with the optimal
approach is insignificant.

##  Appendix A Proof of Theorem 1

In step 2 of the $n$th iteration, we solve Problem (19) to obtain the optimal
$\\{\alpha_{k}^{(n)}\\}_{k\in\cal U}$ and ${\bf{V}}^{(n)}$ with given
${\bf{U}}^{(n-1)}$ and ${\bf{W}}^{(n-1)}$. Hence, we have
${h_{k}}\left({{\bf{V}}^{(n)},{{\bf{U}}_{k}^{(n-1)}},{{\bf{W}}_{k}^{(n-1)}}}\right)\geq\left({\alpha_{k}^{(n)}}\right)^{2}{R_{k,{\rm{min}}}},\forall
k$. In step 3 of the $n$th iteration, we update ${\bf{U}}^{(n)}$ and
${\bf{W}}^{(n)}$ as in (15) with ${\bf{V}}^{(n)}$. According to Lemma 1, we
have
$R_{k}\left({{{\bf{V}}^{(n)}}}\right)={h_{k}}\left({{{\bf{V}}^{(n)}},{\bf{U}}_{k}^{(n)},{\bf{W}}_{k}^{(n)}}\right)\geq{h_{k}}\left({{{\bf{V}}^{(n)}},{\bf{U}}_{k}^{(n-1)},{\bf{W}}_{k}^{(n-1)}}\right)$.
Hence, we have

${h_{k}}\left({{\bf{V}}^{(n)},{{\bf{U}}_{k}^{(n)}},{{\bf{W}}_{k}^{(n)}}}\right)\geq{\left(\alpha_{k}^{(n)}\right)}^{2}{R_{k,{\rm{min}}}}.$
(A.1)

In step 2 of the $(n+1)$th iteration, we obtain
$\\{\alpha_{k}^{(n+1)}\\}_{k\in\cal U}$ and ${\bf{V}}^{(n+1)}$ with given
${\bf{U}}^{(n)}$ and ${\bf{W}}^{(n)}$ by solving Problem (19). Then we have
$\sum\nolimits_{k\in{\cal
U}}{{{\left({\alpha_{k}^{(n+1)}-1}\right)}^{2}}}\leq\sum\nolimits_{k\in{\cal
U}}{{{\left({\alpha_{k}^{(n)}-1}\right)}^{2}}}$. The reason is that from
(A.1), $\\{\alpha_{k}^{(n)}\\}_{k\in\cal U}$ and ${\bf{V}}^{(n)}$ is just a
feasible solution for Problem (19) with given ${\bf{U}}^{(n)}$ and
${\bf{W}}^{(n)}$. Hence, the objective value of Problem (11) is monotonically
decreasing. Obviously, the objective value is lower bounded by zero. Hence,
Algorithm 2 will converge.

##  Appendix B Proof of Theorem 2

We first prove that the sequence of ${\bf{V}}$ generated by the WMMSE
algorithm (i.e. Algorithm 4) always satisfies the rate requirements of Problem
(23). In step 2, we obtain ${\bf{V}}^{(l)}$ with ${\bf{U}}^{(l-1)}$ and
${\bf{W}}^{(l-1)}$. Hence,
${h_{k}}\left({{\bf{V}}^{(l)},{\bf{U}}_{k}^{(l-1)},{\bf{W}}_{k}^{(l-1)}}\right)\geq{R_{k,{\text{min}}}},\forall
k$ hold since ${\bf{V}}^{(l)}$ is feasible for Problem (24). According to
Lemma 1,
${h_{k}}\left({{\bf{V}}^{(l)},{\bf{U}}_{k}^{(l-1)},{\bf{W}}_{k}^{(l-1)}}\right)$
is a lower-bound of ${R_{k}}({\bf{V}}^{(l)})$, i.e.,
${R_{k}}({\bf{V}}^{(l)})\geq{h_{k}}\left({{\bf{V}}^{(l)},{\bf{U}}_{k}^{(l-1)},{\bf{W}}_{k}^{(l-1)}}\right)$.
Hence, ${R_{k}}({\bf{V}}^{(l)})\geq R_{k,{\text{min}}}$ holds. Thus, the
sequence of ${\bf{V}}$ generated by the WMMSE algorithm satisfies the rate
requirements of Problem (23).

Next, we show that the value of the objective function of Problem (23)
monotonically decreases during the iterative process of the WMMSE algorithm.
Denote ${\rm{Obj}}({\bf{V}}^{(l)})$ as the objective value of Problem (23)
when ${\bf{V}}={\bf{V}}^{(l)}$. Step 2 of the WMMSE algorithm updates
${\bf{V}}^{(l)}$ by solving Problem (27) with ${\bf{U}}^{(l-1)}$ and
${\bf{W}}^{(l-1)}$. The objective value of this step,
${\rm{Obj}}({\bf{V}}^{(l)})$, will be no larger than
${\rm{Obj}}({\bf{V}}^{(l-1)})$, i.e.,
${\rm{Obj}}({\bf{V}}^{(l)})\leq{\rm{Obj}}({\bf{V}}^{(l-1)})$. The reason is
that ${\bf{V}}^{(l-1)}$ is a feasible solution for Problem (27) with
${\bf{U}}_{k}^{(l-1)}$ and ${\bf{W}}_{k}^{(l-1)}$ since
${h_{k}}\left({{\bf{V}}^{(l-1)},{\bf{U}}_{k}^{(l-1)},{\bf{W}}_{k}^{(l-1)}}\right)={R_{k}}({\bf{V}}^{(l-1)})\geq{R_{k,{\text{min}}}}$
holds as proved above. In step 3 of the WMMSE algorithm, we update
${\bf{U}}^{(l)}$ and ${\bf{W}}^{(l)}$ by using (15) with ${\bf{V}}^{(l)}$.
This step increases the value of
${h_{k}}\left({{\bf{V}},{\bf{U}}_{k},{\bf{W}}_{k}}\right)$ while maintaining
the same objective value of Problem (23). Therefore, this step provides “room”
for the next iteration to decrease the objective value. In addition, the
objective value is lower bounded by zero. Hence, the WMMSE algorithm
converges.

Then, we prove that given the initial set of precoders, the WMMSE algorithm
converges to a unique solution. Obviously, when ${\bf{V}}$ is given,
${\bf{U}}$ and ${\bf{W}}$ can be uniquely determined by (15). The remaining
task is to prove that given ${\bf{U}}$ and ${\bf{W}}$, the BCD algorithm can
obtain the unique globally optimal solution ${\bf{V}}$. Since
$\\{{\bf{G}}_{k},\forall k\\}$ are positive definite matrices, the objective
function in Problem (27) is a strictly convex function with respect to
(w.r.t.) $\bf{V}$. Obviously, the constraints in Problem (27) are convex
w.r.t. $\bf{V}$ [36]. Hence, Problem (27) is a strictly convex problem [36].
According to [Page 137 in [36]], the globally optimal solution of Problem (27)
is unique. On the other hand, Theorem 3 proves that the BCD algorithm can
obtain the globally optimal solution to the dual problem (32). As Problem (27)
is a convex problem and it satisfies the Slater’s condition [36], the duality
gap between Problem (27) and its dual problem (32) is zero [36]. As a result,
the BCD algorithm can obtain the unique globally optimal solution ${\bf{V}}$.
Finally, by alternatively updating step 2 and step 3, the WMMSE algorithm will
converge to a unique solution. It should be emphasized that as Problem (24) is
non-convex, it may have many locally optimal solutions, and the unique
solution of the WMMSE algorithm depends on the initial point. However, given
the initial points of precoders, the WMMSE algorithm will converge to a unique
solution.

Finally, we prove that the unique solution satisfies the KKT conditions of
Problem (23). Denote the converged solution of the WMMSE algorithm as
${\bf{V}}^{\star}$, ${\bf{U}}^{\star}$ and ${\bf{W}}^{\star}$. With given
${\bf{U}}^{\star}$ and ${\bf{W}}^{\star}$, the Lagrange function of Problem
(24) can be written as

$\displaystyle{\cal L}\left({{\bf{V}},{\bm{\lambda}},{\bm{\mu}}}\right)$
$\displaystyle=$ $\displaystyle\sum\nolimits_{k\in{\cal
U}}{{\bf{\bar{V}}}_{k}^{H}}{{\bf{G}}_{k}}{{{\bf{\bar{V}}}}_{k}}+\sum\nolimits_{k\in\cal
U}{{\lambda_{k}}\left({R_{k,{\rm{min}}}}-{h_{k}}\left({{\bf{V}},{\bf{U}}_{k}^{\star},{\bf{W}}_{k}^{\star}}\right)\right)}$
(B.1) $\displaystyle+$ $\displaystyle\sum\nolimits_{i\in\cal
I}{{\mu_{i}}\left(\sum\nolimits_{k\in{{\cal
U}_{i}}}{\left\|{{{\bf{B}}_{i,k}}{{{\bf{\bar{V}}}}_{k}}}\right\|}_{F}^{2}-{P_{i,\max}}\right)},$

where ${{\bm{\lambda}}}=\left\\{{\lambda_{k},\forall k\in\cal U}\right\\}$ and
${\bm{\mu}}=\left\\{{\mu_{i},\forall i\in\cal I}\right\\}$ are the
corresponding Lagrange multipliers.

According to Theorem 3, the BCD algorithm can obtain the globally optimal
solution of Problem (27) (also Problem (24)) with given ${\bf{U}}^{\star}$ and
${\bf{W}}^{\star}$, there must exist ${{\bm{\lambda}}^{\star}}$ and
${\bm{\mu}^{\star}}$ such that
$\\{{\bf{V}}^{\star},{\bm{\lambda}}^{\star},\bm{\mu}^{\star}\\}$ satisfy the
following KKT conditions

$\displaystyle{\nabla_{{{\bf{\bar{V}}}_{k}}}}{\cal L}$ $\displaystyle=$
$\displaystyle\nabla_{{{\bf{\bar{V}}}_{k}}}\sum\nolimits_{k\in{\cal
U}}{{\bf{\bar{V}}}_{k}^{\star,H}}{{\bf{G}}_{k}}{{{\bf{\bar{V}}}}_{k}^{\star}}-\sum\nolimits_{k\in\cal
U}{\lambda_{k}^{\star}{\nabla_{{{\bf{\bar{V}}}_{k}}}}{h_{k}}\left({{\bf{V}}^{\star},{\bf{U}}_{k}^{\star},{\bf{W}}_{k}^{\star}}\right))}$
(B.6) $\displaystyle+\sum\nolimits_{i\in\cal
I}{\mu_{i}^{\star}}{\nabla_{{{\bf{\bar{V}}}_{k}}}}\left(\sum\nolimits_{k\in{{\cal
U}_{i}}}{\left\|{{{\bf{B}}_{i,k}}{{{\bf{\bar{V}}}}_{k}^{\star}}}\right\|}_{F}^{2}\right)={\bf{0}},\forall
k\in\cal U,$
$\displaystyle\lambda_{k}^{\star}\left({{h_{k}}\left({{\bf{V}}^{\star},{\bf{U}}_{k}^{\star},{\bf{W}}_{k}^{\star}}\right)-{R_{k,{\rm{min}}}}}\right)=0,\forall
k\in\cal U,$
$\displaystyle\mu_{i}^{\star}\left(P_{i,\max}-\sum\nolimits_{k\in{{\cal
U}_{i}}}{\left\|{{{\bf{B}}_{i,k}}{{{\bf{\bar{V}}}}_{k}^{\star}}}\right\|}_{F}^{2}\right)=0,\forall
i\in\cal I,$
$\displaystyle{h_{k}}\left({{\bf{V}}^{\star},{\bf{U}}_{k}^{\star},{\bf{W}}_{k}^{\star}}\right)\geq{R_{k,{\rm{min}}}},\forall
k\in\cal U,$ $\displaystyle\sum\nolimits_{k\in{{\cal
U}_{i}}}{\left\|{{{\bf{B}}_{i,k}}{{{\bf{\bar{V}}}}_{k}^{\star}}}\right\|}_{F}^{2}\leq{P_{i,\max}},\forall
i\in\cal I.$

Since ${\bf{U}}^{\star}$ and ${\bf{W}}^{\star}$ are updated by using (15), we
have
${h_{k}}\left({{\bf{V}}^{\star},{\bf{U}}_{k}^{\star},{\bf{W}}_{k}^{\star}}\right)={R_{k}}({\bf{V}}^{\star})$
according to Lemma 1. By substituting it into the equations (B.6), (B.6) and
(B.6), we find that the set of equations (B.6)-(B.6) are just the KKT
conditions of Problem (23).

##  Appendix C Proof of Theorem 3

According to [36], the dual problem of any optimization problem is a convex
problem. Thus, the dual problem (32) is jointly convex with respect to
${\bm{\lambda}}$ and $\bm{\mu}$. Assuming that the constraint of this problem
satisfies the Slater’s condition, the KKT condition of this problem is
sufficient and necessary for optimality. For given $\bm{\mu}$, the dual
problem (32) is a convex problem w.r.t. ${\bm{\lambda}}$. According to [36],
Newton’s method can obtain the globally optimal solution of dual problem (32)
for given $\bm{\mu}$. In addition, for given ${\bm{\lambda}}$, the dual
problem (32) is convex w.r.t. $\bm{\mu}$, and the gradient descent method can
be applied to obtain the globally optimal solution. Then by adopting the same
idea as in the proof of Theorem 1 in [49], we can prove that the converged
solution also satisfies the KKT condition of Problem (32). Since Problem (32)
is a convex optimization problem, Algorithm 6 can attain the globally optimal
solution of Problem (32).

## References

  * [1] C. Pan, H. Zhu, N. J. Gomes, and J. Wang, “Joint precoding and rrh selection for green mimo c-ran,” in _2016 IEEE Global Communications Conference (GLOBECOM)_ , 2016, pp. 1–5.
  * [2] H. Zhu and J. Wang, “Chunk-based resource allocation in ofdma systems - part i: chunk allocation,” _IEEE Transactions on Communications_ , vol. 57, no. 9, pp. 2734–2744, 2009.
  * [3] ——, “Chunk-based resource allocation in ofdma systems part ii: Joint chunk, power and bit allocation,” _IEEE Transactions on Communications_ , vol. 60, no. 2, pp. 499–509, 2012.
  * [4] H. Zhu, “Radio resource allocation for ofdma systems in high speed environments,” _IEEE Journal on Selected Areas in Communications_ , vol. 30, no. 4, pp. 748–759, 2012.
  * [5] C. Mobile, “C-ran: the road towards green RAN,” _White Paper, ver_ , vol. 2, 2011.
  * [6] N. Alliance, “Suggestions on potential solutions to C-RAN,” _White Paper, January_ , 2013.
  * [7] H. Zhu, “Performance comparison between distributed antenna and microcellular systems,” _IEEE J. Sel. Areas Commun._ , vol. 29, no. 6, pp. 1151–1163, 2011\.
  * [8] J. Wang, H. Zhu, and N. Gomes, “Distributed antenna systems for mobile communications in high speed trains,” _IEEE J. Sel. Areas Commun._ , vol. 30, no. 4, pp. 675–683, 2012.
  * [9] G. Fettweis and E. Zimmermann, “Ict energy consumption-trends and challenges,” in _Proceedings of the 11th International Symposium on Wireless Personal Multimedia Communications_ , vol. 2, no. 4. (Lapland, 2008, p. 6.
  * [10] J. Andrews, S. Buzzi, W. Choi, S. Hanly, A. Lozano, A. Soong, and J. Zhang, “What will 5G be?” _IEEE J. Sel. Areas Commun._ , vol. 32, no. 6, pp. 1065–1082, 2014.
  * [11] C. Pan, W. Xu, J. Wang, H. Ren, W. Zhang, N. Huang, and M. Chen, “Pricing-based distributed energy-efficient beamforming for miso interference channels,” _IEEE Journal on Selected Areas in Communications_ , vol. 34, no. 4, pp. 710–722, 2016.
  * [12] ——, “Totally distributed energy-efficient transmission in mimo interference channels,” _IEEE Transactions on Wireless Communications_ , vol. 14, no. 11, pp. 6325–6338, 2015.
  * [13] C. Pan, W. Xu, W. Zhang, J. Wang, H. Ren, and M. Chen, “Weighted sum energy efficiency maximization in ad hoc networks,” _IEEE Wireless Communications Letters_ , vol. 4, no. 3, pp. 233–236, 2015.
  * [14] H. Ren, N. Liu, C. Pan, and C. He, “Energy efficiency optimization for mimo distributed antenna systems,” _IEEE Transactions on Vehicular Technology_ , vol. PP, no. 99, pp. 1–1, 2016.
  * [15] L. Correia, D. Zeller, O. Blume, D. Ferling, Y. Jading, I. Godor, G. Auer, and L. Van der Perre, “Challenges and enabling technologies for energy aware mobile radio networks,” _IEEE Commun. Mag._ , vol. 48, no. 11, pp. 66–72, 2010.
  * [16] Y. Cheng, M. Pesavento, and A. Philipp, “Joint network optimization and downlink beamforming for CoMP transmissions using mixed integer conic programming,” _IEEE Trans. Signal Process._ , vol. 61, no. 16, pp. 3972–3987, 2013.
  * [17] R. Ramamonjison, A. Haghnegahdar, and V. Bhargava, “Joint optimization of clustering and cooperative beamforming in green cognitive wireless networks,” _IEEE Trans. Wireless Commun._ , vol. 13, no. 2, pp. 982–997, 2014.
  * [18] D. Ng and R. Schober, “Secure and green SWIPT in distributed antenna networks with limited backhaul capacity,” _IEEE Trans. Wireless Commun._ , vol. 14, no. 9, pp. 5082–5097, 2015.
  * [19] F. Zhuang and V. Lau, “Backhaul limited asymmetric cooperation for MIMO cellular networks via semidefinite relaxation,” _IEEE Trans. Signal Process._ , vol. 62, no. 3, pp. 684–693, 2014.
  * [20] V. Ha, L. Le, and N. Dao, “Coordinated multipoint (CoMP) transmission design for Cloud-RANs with limited fronthaul capacity constraints,” _IEEE Trans. Veh. Technol._ , vol. PP, no. 99, pp. 1–1, 2015.
  * [21] M. Tao, E. Chen, H. Zhou, and W. Yu, “Content-centric sparse multicast beamforming for cache-enabled cloud RAN,” _IEEE Trans. Wireless Commun._ , vol. 15, no. 9, pp. 6118–6131, 2016.
  * [22] Y. Shi, J. Zhang, and K. Letaief, “Group sparse beamforming for green Cloud-RAN,” _IEEE Trans. Wireless Commun._ , vol. 13, no. 5, pp. 2809–2823, 2014.
  * [23] S. Luo, R. Zhang, and T. J. Lim, “Downlink and uplink energy minimization through user association and beamforming in C-RAN,” _IEEE Trans. Wireless Commun._ , vol. 14, no. 1, pp. 494–508, 2015.
  * [24] J. Yan, J. Li, L. Zhao, and R. Chen, “Robust joint transmit beamforming with QoS guarantees in time-asynchronous DAS,” _IEEE Trans. Veh. Technol._ , vol. 64, no. 4, pp. 1506–1518, 2015.
  * [25] J. Tang, W. P. Tay, and T. Quek, “Cross-layer resource allocation with elastic service scaling in cloud radio access network,” _IEEE Trans. Wireless Commun._ , vol. 14, no. 9, pp. 5068–5081, 2015.
  * [26] B. Dai and W. Yu, “Energy efficiency of downlink transmission strategies for cloud radio access networks,” _IEEE J. Sel. Areas Commun._ , vol. 34, no. 4, pp. 1037–1050, 2016.
  * [27] E. J. Candes, M. B. Wakin, and S. P. Boyd, “Enhancing sparsity by reweighted $l$1 minimization,” _Journal of Fourier analysis and applications_ , vol. 14, no. 5-6, pp. 877–905, 2008.
  * [28] J. Zhao, T. Quek, and Z. Lei, “Coordinated multipoint transmission with limited backhaul data transfer,” _IEEE Trans. Wireless Commun._ , vol. 12, no. 6, pp. 2762–2775, 2013.
  * [29] T. Huang, Y. Yu, and L. Yi, “Design of highly isolated compact antenna array for MIMO applications,” _International Journal of Antennas and Propagation_ , vol. 2014, 2014.
  * [30] K. Wang, R. Mauermayer, L. Li, and T. Eibert, “A highly compact broadband near-edge antenna for low profile communication devices,” in _Antennas and Propagation (EuCAP), 2015 9th European Conference on_ , 2015, pp. 1–5.
  * [31] Q. Shi, M. Razaviyayn, Z.-Q. Luo, and C. He, “An iteratively weighted MMSE approach to distributed sum-utility maximization for a MIMO interfering broadcast channel,” _IEEE Trans. Signal Process._ , vol. 59, no. 9, pp. 4331–4340, 2011.
  * [32] S. He, Y. Huang, H. Wang, S. Jin, and L. Yang, “Leakage-aware energy-efficient beamforming for heterogeneous multicell multiuser systems,” _IEEE J. Sel. Areas Commun._ , vol. 32, no. 6, pp. 1268–1281, 2014.
  * [33] Y. Li, Y. Tian, and C. Yang, “Energy-efficient coordinated beamforming under minimal data rate constraint of each user,” _IEEE Trans. Veh. Technol._ , vol. 64, no. 6, pp. 2387–2397, 2015.
  * [34] B. Dai and W. Yu, “Backhaul-aware multicell beamforming for downlink cloud radio access network,” in _Communication Workshop (ICCW), 2015 IEEE International Conference on_ , 2015, pp. 2689–2694.
  * [35] ——, “Sparse beamforming and user-centric clustering for downlink cloud radio access network,” _IEEE Access,_ , vol. 2, pp. 1326–1339, 2014.
  * [36] S. Boyd and L. Vandenberghe, _Convex optimization_. Cambridge university press, 2004.
  * [37] W. Nie, F. C. Zheng, X. Wang, W. Zhang, and S. Jin, “User-centric cross-tier base station clustering and cooperation in heterogeneous networks: Rate improvement and energy saving,” _IEEE J. Sel. Areas Commun._ , vol. 34, no. 5, pp. 1192–1206, 2016.
  * [38] W. Feng, N. Ge, and J. Lu, “Hierarchical transmission optimization for massively dense distributed antenna systems,” _IEEE Communications Letters_ , vol. 19, no. 4, pp. 673–676, 2015.
  * [39] T. M. Cover and J. A. Thomas, _Elements of information theory_. John Wiley & Sons, 2012.
  * [40] M. S. Lobo, L. Vandenberghe, S. Boyd, and H. Lebret, “Applications of second-order cone programming,” _Linear algebra and its applications_ , vol. 284, no. 1, pp. 193–228, 1998.
  * [41] B. K. Sriperumbudur, D. A. Torres, and G. R. Lanckriet, “A majorization-minimization approach to the sparse generalized eigenvalue problem,” _Machine learning_ , vol. 85, no. 1-2, pp. 3–39, 2011.
  * [42] D. P. Bertsekas, “Nonlinear programming,” 1999.
  * [43] J. R. Magnus and H. Neudecker, “Matrix differential calculus with applications in statistics and econometrics,” 1995.
  * [44] S. Ye and R. Blum, “Optimized signaling for MIMO interference systems with feedback,” _IEEE Trans. Signal Process._ , vol. 51, no. 11, pp. 2839–2848, 2003.
  * [45] D. Coppersmith and S. Winograd, “Matrix multiplication via arithmetic progressions,” in _Proceedings of the nineteenth annual ACM symposium on Theory of computing_. ACM, 1987, pp. 1–6.
  * [46] L. Venturino, N. Prasad, and X. Wang, “Coordinated linear beamforming in downlink multi-cell wireless networks,” _IEEE Trans. Wireless Commun._ , vol. 9, no. 4, pp. 1451–1461, 2010.
  * [47] G. Auer, V. Giannini, C. Desset, I. Godor, P. Skillermark, M. Olsson, M. Imran, D. Sabella, M. Gonzalez, O. Blume, and A. Fehske, “How much energy is needed to run a wireless network?” _IEEE Wireless Commun. Mag._ , vol. 18, no. 5, pp. 40–49, 2011.
  * [48] A. Dhaini, P.-H. Ho, G. Shen, and B. Shihada, “Energy efficiency in TDMA-based next-generation passive optical access networks,” _IEEE/ACM Transactions on Networking,_ , vol. 22, no. 3, pp. 850–863, 2014\.
  * [49] W. Yu, W. Rhee, S. Boyd, and J. Cioffi, “Iterative water-filling for Gaussian vector multiple-access channels,” _IEEE Trans. Inf. Theory_ , vol. 50, no. 1, pp. 145–152, 2004.