1

Optimization of Massive Full-Dimensional MIMO
for Positioning and Communication

arXiv:1709.02969v1 [] 9 Sep 2017

Seongah Jeong, Osvaldo Simeone, and Joonhyuk Kang

Abstract

Massive Full-Dimensional multiple-input multiple-output (FD-MIMO) base stations (BSs) have the potential to
bring multiplexing and coverage gains by means of three-dimensional (3D) beamforming. Key technical challenges
for their deployment include the presence of limited-resolution front ends and the acquisition of channel state
information (CSI) at the BSs. This paper investigates the use of FD-MIMO BSs to provide simultaneously highrate data communication and mobile 3D positioning in the downlink. The analysis concentrates on the problem of
beamforming design by accounting for imperfect CSI acquisition via Time Division Duplex (TDD)-based training
and for the finite resolution of analog-to-digital converter (ADC) and digital-to-analog converter (DAC) at the BSs.
Both unstructured beamforming and a low-complexity Kronecker beamforming solution are considered, where for
the latter the beamforming vectors are decomposed into separate azimuth and elevation components. The proposed
algorithmic solutions are based on Bussgang theorem, rank-relaxation and successive convex approximation (SCA)
methods. Comprehensive numerical results demonstrate that the proposed schemes can effectively cater to both data
communication and positioning services, providing only minor performance degradations as compared to the more
The work of O. Simeone was partially supported by the U.S. NSF through grant CCF-1525629. O. Simeone has also received funding
from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme (grant agreement
No 725731). This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government(MSIP)
(No. 2017R1A2B2012698).
Seongah Jeong is with the School of Engineering and Applied Sciences (SEAS), Harvard University, Cambridge, MA 02138, USA (Email:
seongah116@gmail.com).
Osvaldo Simeone is with the Department of Informatics, Kings College London, London, UK (Email: osvaldo.simeone@kcl.ac.uk).
Joonhyuk Kang is with the Department of Electrical Engineering, Korea Advanced Institute of Science and Technology (KAIST) Daejeon,
South Korea (Email: jhkang@ee.kaist.ac.kr).

2

conventional cases in which either function is implemented. Moreover, the proposed low-complexity Kronecker
beamforming solutions are seen to guarantee a limited performance loss in the presence of a large number of BS
antennas.

Index Terms

3D beamforming, localization, full-dimensional MIMO (FD-MIMO), digital-to-analog converter (DAC), analogto-digital converter (ADC), Bussgang theorem, successive convex approximation (SCA).

I. I NTRODUCTION
Mobile broadband communication in the New Radio (NR) physical layer of 5G systems is not only
expected to increase data rate and reliability, but also to cater to new services, including proactive radio
resource management, intelligent traffic systems, autonomous vehicles, Internet of Things (IoT), and
device-to-device communication for disaster response and emergency relief. Services such as these can
benefit from location awareness at the mobile users [1], [2]. According to [3], 5G is envisioned to attain
positioning accuracy of one meter or less, outperforming existing positioning techniques such as GPS and
wireless local area network (WLAN) fingerprinting-based systems.
The support for positioning in Long-Term Evolution (LTE) systems has been standardized in the form
of downlink observed time difference of arrival (OTDOA) in Release 9 and unlink TDOA (UTDOA)
in Release 11, with additional work on dedicated signals, procedures, and requirements for vertical
localization accuracy to be carried out for in Release 14, which marks the start of 5G [4].
Among the key technologies introduced to boost the spectral efficiency of 5G, the use of very large
antenna arrays, or massive multiple-input multiple-output (MIMO) systems, [5]–[12] at the base stations
(BSs) is of notable importance. Technical issues that challenge the deployment of massive MIMO systems
include the large space occupation, the hardware cost associated with radio frequency (RF) elements, and
the power dissipation of a large antenna array. As a specific solution, three-dimensional (3D) MIMO,
or Full-Dimensional MIMO (FD-MIMO), has been intensely discussed in the LTE Release 13 [5], [6],

3

[9], [10]. FD-MIMO BSs are equipped with two-dimensional (2D) antenna arrays, thereby reducing the
spatial size of the BS and providing the additional degree of freedom for beamforming design, given by
the elevation angle. In addition, the problem of cost and circuit power dissipation is typically addressed
by using low-resolution Analog-to-Digital converter/Digital-to-Analog converter (ADC/DAC) [7], [8] or
by developing hybrid analog-digital transceiver [13].
In this paper, as illustrated in Fig. 1, we consider a cellular system with FD-MIMO BSs having perantenna limited-resolution front ends. In the system, as per the frame structure in Fig. 2, each MS estimates
its 3D position, as well as the downlink channels, based on the pilot signals received during the downlink
training phase. It then decodes the data received from all the BSs during downlink data phase. It is
noted that downlink localization may be advantageous with respect to uplink localization due to the
large transmission power of the BSs [14]. We focus on the problem of downlink beamforming design
in the presence of imperfect channel state information (CSI), which is estimated at the BSs via uplink
training using time-division duplex (TDD). Unlike the prior work [5], [6], transmit downlink beamforming
is optimized so as to serve both data communication and localization services. The design accounts for
impairments in CSI acquisition accuracy and downlink transmission caused by low-resolution ADC/DACs
at the BSs. Furthermore, we consider both general unstructured beamforming and a low-complexity
Kronecker beamforming solution. In the latter case, the beamforming vectors are decomposed into the
beams in azimuth and elevation and can be optimized separately [12]. The proposed schemes leverage the
Bussgang theorem [15] to model the effects of ADC/DAC quantization. In both cases, we consider two
complementary formulations: (i) sum-rate maximization under localization accuracy and total transmit
power constraints; or (ii) sum-localization error minimization under data rate and total transmit power
constraints. We finally note that our prior work [16] also tackles the problem of beamforming design for
localization and data transmission, but it considers 2D localization and infinite resolution front-ends, and
it concentrates solely on the problem of power minimization.
The rest of the paper is organized as follows. We describe the system model in Section II and

4

𝑁𝐴,𝑗
𝑁𝐸,𝑗
…
…
…
Δz

…
…
…

𝒑𝐵,𝑗

…
…
…

…
…
…
MS i
𝒑𝑀,𝑖 =[𝑥𝑀,𝑖 𝑦𝑀,𝑖 𝑧𝑀,𝑖 ]𝑇

𝜙𝑗𝑖

BS j
=[𝑥𝐵,𝑗 𝑦𝐵,𝑗 𝑧𝐵,𝑗 ]𝑇
…
…
…

𝜃𝑗𝑖

Δx

𝑑𝑗𝑖

z
y
x

Δy

Fig. 1. In the system under study, NB BSs equipped with FD-MIMO arrays and a limited-resolution front end serve NM single-antenna
MSs. Each MS estimates its 3D position p M,i based on the pilot signals received during the downlink training phase and receives data from
all the BSs during the downlink data phase. The system operates via TDD and uplink and downlink channels are estimated at the BSs during
the uplink training phase.

the performance metrics in Section III. Then, we propose unstructured beamforming and Kronecker
beamforming designs in Section IV and Section V, respectively. In Section VI, numerical results are
presented, and concluding remarks are summarized in Section VII.

II. S YSTEM M ODEL
Consider the cellular system as illustrated in Fig. 1, in which NB BSs, equipped with massive FDMIMO antenna arrays and a limited-resolution front-end, serve NM MSs for the purpose of downlink
communication and positioning. As detailed below, the system operates via TDD for the purpose of
uplink-based channel estimation at the BSs [17], and it uses a frequency reuse scheme that assigns a
different band to each BS in the area under study (see, e.g,. [18]). Our focus is the design of downlink
beamforming vectors based on estimated CSI with the goal of ensuring performance guarantees in terms
of both data transmission and localization in the presence of finite-resolution front-ends at the BSs.
The set of BSs and MSs are denoted as NB = {1, . . . , NB } and NM = {1, . . . , NM }, respectively. The
FD, or 2D, uniform rectangular array (URA) at each BS j ∈ NB has Nj = NA,j NE,j antennas, where
the NA,j horizontal antennas, placed along the y axis, have spacing dA , and the NE,j vertical antennas,

5

Coherence time

①
Uplink pilot
(𝑛ത𝑝 symbols)
Channel estimation
at the CU

Uplink data

②
Downlink pilot
(𝑛𝑝 symbols)
Localization
& channel estimation
at the MSs

③

Downlink data
(𝑛𝑑 symbols)
Data decoding
at the MSs

Fig. 2. Illustration of the frame structure in each frequency band.

located along the z axis, have spacing dE . Note that, as in [11], [12], we assume no mechanical downtilt
for the antenna array. The MSs have a single antenna.
We assume a digital massive MIMO implementation in which a low-resolution ADC/DAC is available
for each antenna element at the BSs [7], [8], [19]–[21]. Each ADC/DAC at BS j has Bj quantization levels
Bj = {bj,0 , . . . , bj,Bj −1 } for both the in-phase and quadrature components. We define the corresponding
quantization function operating separately on the in-phase and quadrature components of each element of
the argument vector as QBj (·).
The MS i ∈ NM is located at position p M,i = [xM,i yM,i zM,i ]T , which is randomly and uniformly
distributed within a ∆x × ∆y × ∆z cube. Instead, BS j ∈ NB is located at a fixed position p B,j =
[xB,j yB,j zB,j ]T within the cube. The positions {ppB,j }j∈NB of the BSs are assumed to be known to all the
nodes in the network. The distance dji , azimuth angle φji and elevation angle θji between BS j and MS i are

 q
2
2


(xM,i −xB,j ) +(yM,i −yB,j )
x
−x
M,i
B,j
−1
−1
− yM,i −yB,j , and θji = tan
denoted as dji = kppM,i −ppB,j k, φji = tan
,
−
zM,i −zB,j
respectively, where the angles φji and θji are defined with respect to the negative y-axis and z-axis.
As in [16], we assume that each BS j communicates with the MSs simultaneously over a frequency
band j that is orthogonal to the bands i assigned to any other BS i 6= j. Note that this precludes the use of
cooperative processing across BSs, such as cooperative multi-point (CoMP) transmission [18]. As shown
in Fig. 2, transmission in each frequency band is organized in frames. Uplink and downlink channels in all
bands are constant within each frame and change according to stationary independent ergodic processes
from one frame to the next. Due to channel reciprocity, the uplink channel matrix is assumed to be equal

6

to the transpose of the downlink channel matrix [17].
Frames are divided into the following slots. 1) Uplink training: The MSs transmit simultaneously
orthogonal pilot signals in all the bands. Each BS estimates the channels of all the MSs in the assigned
band. This pilot slot, of duration n̄p symbols, is followed by an uplink data slot, which is not of interest in
this work. 2) Downlink training: Each MS uses the signal received in the downlink pilot slot to perform
CSI estimation and localization. To enable 3D localization, we assume the condition NB ≥ 4, so that each
MS i can locate the position pM,i via multiangulation or multilateration based on the time measurements.
Furthermore, the BSs and MSs are assumed to have a common time reference. 3) Downlink data: The
MSs receive data in the downlink data slot from all the BSs.
In the following, we will use subscripts or superscripts p and d for variables related to pilots and data,
respectively, while the bar notation, plain letters and hat notation are used for variables pertaining to the
uplink transmission, downlink transmission and estimation, respectively.

A. FD Channel Models
We will consider two types of FD channel models that are typically adopted for 3D modeling, namely
the correlated channel model [9], [10] and the Kronecker channel model [11], [12].
1) Correlated Channel Model: The channel between BS j and MS i in any given coherence time is
denoted as
g ji = ζjih ji ,

(1)

where ζji accounts for the path loss between BS j and MS i, which can be written as

ζji =


1+

dji
d0

η − 21
,

(2)

with η being the path loss exponent and d0 being a reference distance (see, e.g., [11], [16], [22]). The
pathloss (2) is random, given that the MSs are spatially and randomly distributed within the given cube
volume, with first and second moments E[ζji ] = µζj and E[ζji2 ] = σζ2j , respectively. Each vector h ji ∼
CN (00Nj ×1 , R hji ) represents correlated small-scale Rayleigh fading and is characterized by a covariance

7

hjih H
matrix R hji , E[h
ji ]. The path loss and small-scale fading parameters are assumed to be independent,
2
which yields the moments E[gg ji ] = 0 Nj ×1 and E[gg jig H
ji ] = σζj R hji , R gji . Following [9]–[12], the

correlation between the (k, l)th antenna element and (p, q)th antenna element, with the (k, l)th antenna
element indicating the kth in elevation and lth in azimuth element of the URA, is given as
2


γ1 − γ7 j γ2 γ6 − (γ2 σφ2γsin φji )
5
,
R hji (k,l),(p,q) = √ e 2γ5 e γ5 e
γ5

(3)

where {γk }k∈{1,...,7} are functions of (k, l, p, q) and are calculated as in [12] (see Appendix A).
2) Kronecker Channel Model: According to the simplified Kronecker model [11], [12], the channel
vector g ji can be decomposed as
g ji = ζjih A,ji ⊗ h E,ji ,

(4)

where ζji is defined as in (2), while h A,ji ∼ CN (00NA,j ×1 , R hA,ji ) and h E,ji ∼ CN (00NE,j ×1 , R hE,ji ) are the
NA,j × 1 azimuth and NE,j × 1 elevation components, respectively, with covariance matrices R hA,ji ,


hE,jihH
hA,jihH
E[h
E,ji ]. These can be computed as R hA,ji l,q =
A,ji ] and R hE,ji , E[h
−

e

(γ2 σφ sin φji )2
2γ5

−
√1 e
γ5

2 cos2 φ
γ3
ji
2γ5

e

j

γ2 cos φji
γ5



2πdE 2


2πdE
(p−k)2 sin2 θji
−1 σ
, where [·]l,q denotes the (l, q)th
and R hE,ji k,p = ej λ (p−k) cos θji e 2 θ λ

entry of the argument matrix and the parameters {γk }k∈{1,...,7} are defined in Appendix A [12]. The
Kronecker model has been reported to be a valid approximation for propagation scenarios characterized
by scatters distributed in the vicinity of BS and MS [23]. We note that the correlated channel model
includes the Kronecker channel model as a special case.

B. Uplink Signal Model
In this section, we model the pilot signals received at the BS j ∈ NB during uplink training and the
estimated channels of all the MSs at the BSs by accounting for the low-resolution ADCs of the BSs as
in [7], [19], [20].
In each uplink pilot slot of any band j, the MSs simultaneously transmit orthogonal pilot signals
s̄spji ∈ Cn̄p ×1 = [s̄pji (1) · · · s̄pji (n̄p )]T of duration n̄p ≥ NM symbols with normalized energy ks̄spji k2 /n̄p = 1
for all i ∈ NM .

8

To analyze the impact of quantization resulting from the low-resolution ADCs, as in [7], [19], [20], we
leverage the Bussgang theorem [15]. Accordingly, we write the Nj × n̄p output of the ADCs at BS j as
!
X p
X p
p
p
p
pH
Z j = (1 − Dj )
Qj ,
U j = QBj
P̄gg jis̄sji + Z̄
P̄gg jis̄spH
(5)
Ū
ji + Q̃
i∈NM

i∈NM
p

Z j = [z̄z pj (1) · · · z̄z pj (n̄p )] is the Nj × n̄p noise
where P̄ is a transmit power per MS for uplink training; Z̄
matrix, which consists of independent and identically distributed (i.i.d.) CN (0, N0 ) entries; Dj is tabulated
p

Qj = [q̃q pj (1) · · · q̃q pj (n̄p )] is the Nj × n̄p noise matrix
in [24] for different values of the resolution Bj ; and Q̃
including the channel and quantization noise. Each Nj × 1 noise vector q̃q pj (l) has a diagonal correlation
P
matrix R q̃jp = (1 − Dj )(N0I Nj ×Nj + Dj P̄ diag{ i∈NM R gji }), with diag{·} being a diagonal matrix whose
diagonal elements are the same as the argument matrix, while off-diagonal elements are equal to zero.
p

Qj is uncorrelated but not independent of the desired signal in (5), and it is
The effective noise matrix Q̃
also not Gaussian.
Using the linear minimum mean-square-error (LMMSE) approach, the channel g ji between BS j and
MS i can be estimated based on the received signal (5), yielding the estimate [25]

−1
1
p
√ R gji R gji +
U j s̄spji .
Ū
R q̃jp
ĝg ji =
2
(1 − Dj ) n̄p P̄
(1 − Dj )n̄p P̄
1

(6)

From the orthogonality property of the LMMSE [25], the channel vector g ji can be written as a function of
∆gji , where ĝg ji and ∆ gji are zero-mean and uncorrelated with covariance mathe estimate (6) as g ji = ĝg ji +∆

−1

−1
1
1
p
trices R ĝji = R gji R gji + (1−Dj1)2 n̄p P̄ R q̃jp
RH
and
R
R
R
+
R
RH
.
=
∆
g
g
2
2
gji
gji
ji
ji
q̃ p
(1−Dj ) n̄p P̄
(1−Dj ) n̄p P̄ q̃j
j

C. Downlink Signal Model
Here, we describe the signal model for downlink pilot and data transmission. In a manner similar to
the treatment of uplink training, we take into account the impact of the limited resolution DACs available
at the BSs.

9

In both downlink training and data phases, due to the per-antenna low-resolution DAC at the BS j, the
transmitted signals are quantized before transmission. The discrete-time transmitted signals are and given
as [15]
W j s νj (l)) = (1 − Dj )
u νj (l) = QBj (W

X

w ji sνji (l) + q νj (l),

(7)

i∈NM

where the superscript ν ∈ {p, d} indicates the training, or pilot (p), phase and the data (d) phase; W j =
w j1 · · · w jNM ] ∈ CNj ×NM is the beamforming matrix, with w ji being the Nj × 1 beamforming vector
[w
employed by BS j to communicate with MS i in the allocated band; and sνj (l) = [sνj1 (l) · · · sνjNM (l)]T ∈
CNM ×1 with spji (l), for l = 1, . . . , np , and sdji (l), for l = 1, . . . , nd , are pilot symbols and data symbols,
respectively, which are assumed to be i.i.d. CN (0, 1) variables [26]. In (7), Bussgang theorem [15] is
applied to model quantization, where Dj is defined as in (5), and q νj (l) is the uncorrelated zero-mean
quantization error vector with diagonal correlation matrix
W j ) = Dj (1 − Dj )diag{W
W jW H
R qj (W
j }.

(8)

w ji } are jointly designed based on the estimated CSI vectors
In this work, the beamforming vectors {w
{ĝg ji } in (6). In practice, this design can be carried out at a central unit (not shown in Fig. 1) that collects
w ji },
all the estimates from the BSs. We consider two beamforming design approaches for the vectors {w
w ji };
namely: (i) unstructured design, in which no constraints are imposed on the structure of the vectors {w
(ii) structured design based on the Kronecker parameterization
w ji = w A,ji ⊗ w E,ji ,

(9)

with w A,ji being its NA,j × 1 azimuth component and w E,ji being its NE,j × 1 elevation component. This
parameterization is particularly well studied for the scenarios in which the channel is well described by
the Kronecker channel model and that it has been previously considered in [11], [12]. It has the advantage
that it decreases the number of parameters to be designed from NM
to NM

P

j∈NB (NA,j

P

j∈NB Nj = NM

P

j∈NB

NA,j NE,j

+ NE,j ), hence significantly reducing the complexity of the design problem, specially

in multiuser massive MIMO systems. We will study unstructured beamforming design under the general

10

correlated channel model, and Kronecker beamforming design assuming a Kronecker channel model in
order to obtain lower-complexity solutions.
The BSs’ antennas apply the analog filter p(t) prior to transmission, where p(t) is a unitary-energy
ν
Nyquist pulse whose Fourier transform is P (f ). As a result, the continuous-time signal yji
(t) received at

MS i during the pilot phase for ν = p or during the data phase for ν = d can be expressed as
ν
ν
yji
(t) = g Tjiu νj (t − τji ) + zji
(t),

where uνj (t) =

Pnν

ν
l=1 u j (l)p(t

(10)

− (l − 1)Ts ) is the modulated transmitted signal, with Ts being the symbol

period; τji = dji /c is the effective propagation delay between BS j and MS i, with c being the propagation
ν
speed; and the noise zji
(t) is complex white Gaussian with zero mean and two-sided power spectral density
p
(t)
N0 . Each MS i estimates its position p M,i as well as the downlink channels based on the signals yji
d
received during downlink training phase, and then decodes the signals yji
(t) received during downlink data

phase based on the available CSI obtained via downlink training. The consideration of the continuous-time
signal (10) will be useful below when considering the performance of localization.
ν
Each MS i passes the received signal yji
(t) through a filter matched to p(t). Assuming time synchro-

nization at the MS [26], we can write the discrete-time received signal at MS i in terms of the effective
complex gain between BS j and MS i expressed in terms of the BSs’ channel estimate ĝg ji , which is
defined as

(k)
w jk ) = (1 − Dj )ĝg Tjiw jk .
αji (w

(11)

This is done as follows:
ν
ν
yji
(l) = g Tjiu νj (l) + zji
(l)

= (1 − Dj )

X

ĝg Tjiw jk sνjk (l) + (1 − Dj )

k∈NM

=

X

(k)
ν
w jk )sνjk (l) + z̃ji
αji (w
(l),

(12a)
X

ν
∆ Tgji w jk sνjk (l) + g Tjiq νj (l) + zji
(l)

(12b)

k∈NM

(12c)

k∈NM
ν
ν
for l = 1, . . . , nν , where zji
(l) ∼ CN (0, N0 ) is i.i.d. additive Gaussian noise; and z̃ji
(l) is a non-Gaussian

effective noise term, which includes channel estimation noise, distortion noise and channel noise with

11

power
X

W j) =
σz̃2ji (W

R gji diag{W
W jW H
R ∆gji w jkw H
(1 − Dj )2 tr{R
jk } + Dj (1 − Dj )tr{R
j }} + N0 .

(13)

k∈NM

Note that (12) is written as a function of the CSI available at the BSs, so as to obtain expressions for
the performance metrics of interest (see next section) that can be optimized at a central unit connected
all BSs.

III. P ERFORMANCE M ETRICS
As discussed, we are interested in designing beamforming vectors based on the estimated CSI available
at the BSs so as to guarantee performance requirements in terms of both data transmission rate and
localization accuracy for the downlink. In this section, we discuss the calculation of the performance
criteria of achievable transmission rates and of localization accuracy.

A. Transmission Rate
Each MS i decodes on any band j the signal received from BS j. In order to evaluate achievable
rates that can be used at the BSs for beamforming design, we treat the additive noise and interference
in (12) as Gaussian and independent of the signal. As proved in [27], this yields a lower bound on the
achievable rate based on signal (12). We note that the resulting rate is, strictly speaking, only achievable
(k)
w jk )}i,k∈NM in (12). An additional noise
if the MS can estimate correctly the effective channel gain {αji (w

term could be added in order to account for the channel estimation errors at the MSs, but this is not done
here so as to avoid introducing more notation.
Treating the effective noise term

P

k∈NM ,k6=,i

(k)
d
w jk )sdjk (l)+ z̃ji
αji (w
(l) in (12), which includes interference

d
from the undesired signals as well as the noise z̃ji
(l), as Gaussian and independent of the useful signal,

we obtain the rate

W j) ,
Cji (W
where n = np + nd .

(i)

nd

log2 1 + P
NB n

2



2


 (bps/Hz) ,

w ji )
αji (w

(k)
w jk )
k∈NM ,k6=i αji (w

+

W j)
σz̃2ji (W

(14)

12

B. Localization Accuracy
MS i estimates its position p M,i based on the received pilot signals y pi = [(yy p1i )T · · · , (yy pNB i )T ]T , where
p
y pji is the vector representation of the time-series signal yji
(t) in (10). In order to evaluate the localization

accuracy, we adopt the standard squared position error (SPE) criterion [16], [28]–[30], which is defined
for each MS i as
W , pM,i , ĝg i ) = E
SPE(W

h

p̂pM,i − pM,i

2

i
pM,i , ĝg i ,

(15)

with p̂pM,i being the position estimate at MS i and ĝg i = {ĝg ji }j∈NB being the channel estimates available
at the BSs. The expectation in (15) is conditioned on the unknown parameters including user’s position
and the channel estimates. We evaluate a lower bound on (15) by considering the modified Cramér-Rao
bound (MCRB) [31], [32] that is obtained as the trace of the inverse of the average equivalent Fisher
information matrix (EFIM), when the average is over the pilot sequence. We specifically use the lower
bound of the EFIM derived in [33], [34] by treating the additive noise in (12) as Gaussian. Accordingly,
the resulting performance metric for the SPE of MS i as the MCRB is given as

W , p M,i , ĝg i ) ≤ ρ(W
W , p M,i , ĝg i ) , tr J −1 (W
W , p M,i , ĝg i ) ,
SPE(W
where
2

2

8π np β X
W , p M,i , ĝg i ) =
J (W
c2
j∈N

P

2

8π np β X
=
c2
j∈N

2

j jij Tji

W j)
σz̃2ji (W

B

2

k∈NM

(k)
w jk )
αji (w

P

B

k∈NM

(k)
w jk )
αji (w

W j)
σz̃2ji (W

(16)

2

J φ,θ (φji , θji ),

(17)

R

with β = { |f P (f )|2 df }1/2 being the effective bandwidth and P (f ) being the Fourier transform of the
filter p(t), j ji =

1
(ppM,i
dji

− p B,j ) = [sin φji sin θji cos φji sin θji cos θji ]T and J φ,θ (φ, θ) = j jij Tji . The

detailed derivation of (16) is provided in Appendix B.

C. Total Transmit Power
We conclude this section by evaluating the transmitted power for the BSs as a function of the beamforming vectors. Applying the Bussgang theorem [15] to the quantized precoded signals in (7), the total

13

transmit powers during training phase and data phase are given as
X

W j) =
Pj (W

j∈NB

X

h
E

u νj (l)

2

i

=

X


W jW H
tr (1 − Dj )W
.
j

(18)

j∈NB

j∈NB

Note that the transmit powers at BS j during training phase and data phase are equal.
IV. U NSTRUCTURED B EAMFORMING D ESIGN
In this section, we aim at optimizing the beamforming strategy by assuming the correlated channel
model described in Section II-A1 and without imposing any structure on the beamforming matrix. We
will focus on two dual problems: (i) sum-rate maximization under localization accuracy and total transmit
power constraints; and (ii) sum-SPE minimization under the data rate and total transmit power constraints.
We here consider the average localization performance with respect to the MSs’ positions, given that the
MSs are spatially and randomly distributed within the given cube.
A. Problem Formulations
,
We denote the data rate and localization accuracy requirements for each MS i as Cimin and ρmax
i
respectively, while the constraint of total transmit power expenditure in the downlink is denoted as P .
W j ) in (14), the
By using the performance metrics discussed in Section III, namely the rate function Cji (W
W , p M,i , ĝg i ) in (16) and the transmit power function Pj (W
W j ) in (18), two design problems
SPE function ρ(W
under study are formulated as
(P1: Sum-rate Maximization)

X

maximize
W

W j)
Cji (W

(19a)

j∈NB ,i∈NM

s.t.

W , p M,i , ĝg i )] ≤ ρmax
Ep M,i [ρ(W
, for all i ∈ NM
i

(19b)

X

(19c)

W j ) ≤ P,
Pj (W

j∈NB

and
(P2: Sum-SPE Minimization)

X

minimize
W

s.t.

W , p M,i , ĝg i )]
Ep M,i [ρ(W

(20a)

i∈NM

X

W j ) ≥ Cimin , for all i ∈ NM
Cji (W

(20b)

W j ) ≤ P.
Pj (W

(20c)

j∈NB

X
j∈NB

14

The problems P1 and P2 are challenging to solve due to the expectation in the localization error performance (19b) and (20a) with respect to the unknown users’ locations, and to the non-convexity of both
rate and localization criteria (14) and (16). To address the first issue, we approximate the expectation in
(19b) and (20a) by using the sample average approximation method [35]. Accordingly, we estimate the
W , p M,i , ĝg i )] with respect to the user position as
expectation of Ep M,i [ρ(W
Ns
1 X
W , p M,i,m , ĝg i ),
W , p M,i , ĝg i )] ≈
ρ(W
Ep M,i [ρ(W
Ns m=1

(21)

for i ∈ NM , where pM,i,1 , pM,i,2 , . . . , pM,i,Ns are Ns independent realizations of the MS i’s position. The
latter is assumed to be uniformly distributed in the ∆x × ∆y × ∆z cube as described in Section II 1 . To
tackle the second issue, we apply rank relaxation and the successive convex approximation (SCA) method
introduced in [36], [37] as detailed in the following.

B. Sum-Rate Maximization
In this section, we elaborate on the solution of sum-rate maximization problem P1 with the stochastic
approximation (21). To this end, we start by reformulating the problem with respect to the beamforming
Ωji }j∈NB ,i∈NM with Ω ji = w jiw H
covariance matrices Ω = {Ω
ji . Accordingly, the elements of the diagonal
W j ) of the distortion noise in (8) during both training phase and data phase can
covariance matrix R qj (W
be expressed in terms of Ω as
X


W j ) n,n = Dj (1 − Dj )
E nΩ jiE Tn },
R qj (W
tr{E

(22)

i∈NM

for n = 1, . . . , Nj , where we have defined the matrix E n = eneTn , with en being a vector whose nth entry
(k)

w jk )|2 in (11) can be
equals to 1 and the rest equal to zero. Similarly, the effective power gain |αji,k (w
(k)
Ωjk ) = (1 − Dj )2ĝg H
g ji . Introducing the auxiliary variables κ = {κjn }j∈NB ,n=1,...,Nj for
written as ξji (Ω
ji Ω jkĝ
1

In principle, based on the channel estimates {ĝg i }, the optimizer could restrict the uncertainty area for the user to a smaller volume. The

assumed uniform distribution can hence be thought of as providing a worst-case performance.

15

each BS antenna, where κjn ∈ R and variables χ = {χji }j∈NB ,i∈NM for each BS-MS pair, where χji ∈ R,
a rank-relaxed version of problem P1 based on the empirical approximation (21) can be written as
(P1-1): maximize
κ,χ
χ
Ω ,κ

X
nd
log2
NB n j∈N ,i∈N
B

M

+

X

(k)
Ωjk ) + N0 +
ξji (Ω

k∈NM
Nj
X

X

R ∆gji Ω jk }
(1 − Dj )2 tr{R

k∈NM



n=1

+

X


κjn tr E nR gji E Tn  − log2

X
k∈NM

(k)
Ωjk ) + N0
ξji (Ω

k∈NM ,k6=,i

R ∆gji Ω jk } +
(1 − Dj )2 tr{R

Nj
X



κjn tr E nR gji E Tn 

(23a)

n=1


!−1 
NB
Ns


2
2 X
X
8π np β
1
s.t.
tr
J
(φ
,
θ
)χ
, i ∈ NM
≤ ρmax
φ,θ
ji,m ji,m
ji
i

Ns m=1 
c2
j=1
X
Ωji } ≤ P
(1 − Dj )tr {Ω

(23b)
(23c)

j∈NB ,i∈NM
(k)
Ωjk )
ξji (Ω
≥ χji , j ∈ NB , i ∈ NM

PNj
P
R ∆gji Ω jk } + n=1
κjn tr E nR gji E Tn
N0 + k∈NM (1 − Dj )2 tr{R
X 
Dj (1 − Dj )
tr E nΩ jiE Tn ≤ κjn , j ∈ NB and n = 1, . . . , Nj

P

k∈NM

(23d)
(23e)

i∈NM

Ω ji  0, j ∈ NB , i ∈ NM

(23f)

χji ≥ 0, j ∈ NB , i ∈ NM ,

(23g)

where φji,m and θji,m are azimuth and elevation angle between the Ns independent realizations of the MS
Ωji ) = 1
i’s position and BS j. The rank-relaxation in (23) is obtained by dropping the constraint rank(Ω
for j ∈ NB and i ∈ NM . Note that the equivalence between problem P1 and P1-1 under the rank
relaxation follows from the fact that, at an optimum, inequalities (23d) and (23e) can be assumed to hold
with equality without loss of optimality. The rank-relaxed problem P1-1 is still not convex owing to the
presence of the non-convex objective function (23a) and non-convex constraints (23d). To resolve this
problem, we apply the SCA method introduced in [36], [37], which yields an iterative algorithm that is
guaranteed to converge to a stationary point of the original non-convex problem under suitable conditions.
This algorithm solves a sequence of strongly convex problems obtained as a local approximation of the
original non-convex problem. In order to develop the SCA-based algorithm, we use the following lemmas.

16

x) =
Lemma 1 ([43, Example 7]): Consider a non-concave utility function U (x

PI

i=1

x) with fi (x
x) =
fi (x

x) − fi− (x
x) being a difference of concave (DC) function, where fi+ (x
x) and fi− (x
x) are concave and
fi+ (x
x), a concave approximant of U (x
x) that
continuously differentiable. Then, for any y in the domain of U (x
has the properties required by the SCA algorithm [36, Assumption 2] is given as
x; y ) ,
Û (x

I
X

x) − ∇x fi− (yy )T (x
x − y) −
fi+ (x

i=1

τi
xi − y i )T H i (yy )(x
xi − y i ),
(x
2

(24)

where τi > 0 is a positive constant ensuring that (24) is strongly concave and H i (yy ) is any uniformly
positive definite matrix.
x1 , x2 ) ≤ 0, where g(x
x1 , x2 ) =
Lemma 2 ([43, Example 4]): Consider a non-convex constraint g(x
x1 )h2 (x
x2 ) and h1 (x
x1 ) and h2 (x
x2 ) are convex and non-negative. Then, for any (yy 1 , y 2 ) in the domain
h1 (x
x1 , x 2 ), a convex approximant that satisfies the conditions [36, Assumption 3] required by the SCA
of g(x
algorithm is given as

1
1 2
x1 ) + h2 (x
x2 ))2 −
(h1 (x
h1 (yy 1 ) + h22 (yy 2 )
2
2

x1 , x 2 ; y 1 , y 2 ) ,
ḡ(x

x1 − y 1 ) − h2 (yy 2 )∇x2 hT2 (yy 2 )(x
x2 − y 2 ).
−h1 (yy 1 )∇x1 hT1 (yy 1 )(x

(25)

Using Lemma 1 for the objective function (23a), Appendix C shows that a strongly concave surrogate
function can be obtained as
nd
Ĉ(vv ; v (t)) ,
NB n

!
X

fji+ (vv ) −

X

j,i

∇v fji− (vv (t))T (vv − v (t))

j,i

−

X τΩji
2

j,i

Ωji − Ω ji (t)k2F −
kΩ

X τκjn
j,n

2

(κjn − κjn (t))2 ,

(26)

where fji+ (vv ) and fji− (vv ) are defined as (33);
∇v fji− (vv (t))T (vv − v (t)) =

+

1
ln 2

X

Ωjk − Ω jk (t))ĝg ji
(1 − Dj )2ĝg H
ji (Ω

k∈NM ,k6=,i


Nj
X

Ωjk − Ωjk (t))} +
R ∆gji (Ω
(1 − Dj )2 tr{R
(κjn − κjn (t))tr E nR gji E Tn 

X

n=1

k∈NM


÷

X
k∈NM ,k6=,i

(k)
Ωjk (t)) + N0 +
ξji (Ω

X
k∈NM

R ∆gji Ω jk (t)} +
(1 − Dj )2 tr{R

Nj
X
n=1



κjn (t)tr E nR gji E Tn  ; (27)

17

Ω, κ , χ ) and v (t) = (Ω
Ω(t), κ (t), χ (t)) as the tth
and we have defined the set of optimization variables v = (Ω
iterate of the SCA algorithm. Furthermore, applying Lemma 2 to constraint (23d), we obtain the following
strongly concave approximation of the problem P1-1 for a given a feasible solution v (t) as
(P1-2): maximize
v

Ĉ(vv ; v (t))

(28a)
X

s.t. ḡji (v1 , v 2 ; v1 (t), v 2 (t)) ≤

(k)
Ωjk ) − χji N0 , j ∈ NB , i ∈ NM
ξji (Ω

(28b)

k∈NM

(23b), (23c), (23e)-(23g),

(28c)

where
ḡji (v1 , v 2 ; v1 (t), v 2 (t)) =

1
1 2
(h1,ji (v1 ) + h2,ji (vv 2 ))2 −
h1,ji (v1 (t)) + h22,ji (vv 2 (t)) − h1,ji (v1 (t))(χji − χji (t))
2
2


Nj
X
X

R ∆gji (Ω
Ωjk − Ω jk (t))} +
(κjn − κjn (t))tr E nR gji E Tn  ;
(1 − Dj )2 tr{R
−h2,ji (vv 2 (t)) 
n=1

k∈NM

h1,ji (v1 ) = χji ; h2,ji (vv 2 ) =

(29)

P

k∈NM (1

R ∆gji Ω jk } +
− Dj )2 tr{R

PNj

n=1


κjn tr E nR gji E Tn ; v1 = χji ; v 2 =

Ωj1 , . . . , ΩjNM , κj1 , . . . , κjNj }; v1 (t) = χji (t) and v 2 = {Ω
Ωj1 (t), . . . , ΩjNM (t), κj1 (t), . . . , κjNj (t)}. Ad{Ω
ditional details for the derivation of (28) can be found in Appendix C. The problem P1-2 has a unique
solution denoted by v̂v (vv (t)). Using problem P1-2, the proposed SCA-based algorithm is summarized in
Ωopt , κ opt , χ opt ) using SCA, and then computes a feasible
Algorithm 1. Algorithm 1 obtains a solution (Ω
w opt
Ωopt
beamforming vectors {w
ji } from the covariance matrices {Ω
ji } using standard rank-reduction method
coupled with the scaling method [38]. It is noted that the rate and the EFIM functions are monotonic
with respect to the scaling factor s , which entails that Algorithm 1 provides a feasible solution for the
original problem P1.

C. Sum-SPE Minimization
The sum-SPE minimization problem P2 with the stochastic approximation (21) in lieu of (20a) can be
addressed by using rank relation and SCA in a manner similar to that detailed above for problem P1.
Specifically, by introducing the beamforming covariance matrices Ω and the auxiliary variables κ and χ

18

Algorithm 1 Unstructured beamforming design for sum-rate maximization & sum-SPE minimization
Ω(0), κ (0), χ (0)} and fix parameters τΩji , τκjn > 0 for j ∈ NB ,
Initialization: Initialize v (0) = {Ω
i ∈ NM , n = {1, . . . , Nj } and v , s > 0. Set t = 0.
Repeat
1. Compute v̂v (vv (t)) using P1-2 for sum-rate maximization (or P2-2 for sum-SPE minimization);
2. Set v (t + 1) = v (t) + (t)(v̂v (vv (t)) − v (t)) for some (t) ∈ (0, 1];
3. t ← t + 1;
4. If kvv (t) − v (t − 1)k < v , stop. Otherwise, go to step 1.
end
Ωopt , κ opt , χ opt ) ← v̂v (vv (t)).
(Ω
q
Ωopt
v max (Ω
Ωopt
Ωopt
w ji = λmax (Ω
Rank reduction: Extract ŵ
ji )v
ji ) for all j ∈ NB and i ∈ NM , where λmax (Ω
ji )
opt
Ωji ) are the maximum eigenvalue and the corresponding eigenvector of the beamforming
and v max (Ω
opt
Ω
matrix ji , respectively.
w ji is feasible for the sum-rate maximization (or the sum-SPE minimization)
Scaling: Check whether ŵ
opt
w
w
w ji ← (1 + s )ŵ
w ji for any positive integer s until ŵ
w ji is
or not. If so, ji = ŵ ji . Otherwise, rescale ŵ
feasible.
as in Section IV-B and following similar steps, we derive in Appendix D the SCA algorithm detailed in
Algorithm 1, where the strongly convex problem P2-2 is defined as

!−1 
NB
Ns
 8π 2 n β 2 X

X 1 X
p
tr
(P2-2): minimize
J
(φ
,
θ
)χ
φ,θ
ji,m
ji,m
ji
v

Ns m=1 
c2
j=1
i∈NM
X
s.t.
C̄ji (vv ; v (t)) ≥ Cimin , i ∈ NM

(30a)
(30b)

j∈NB

(28b), (23c), (23e)-(23g).

(30c)

Problem P2-2 has a unique solution denoted by v̂v (vv (t)). In P2-2, we have defined the set of optimization
Ω, κ , χ ) and
variables for the rank-relaxed problem P2-1 in (35) for sum-SPE minimization P2 as v = (Ω
Ω(t), κ (t), χ (t)) for the tth iterate within the feasible set of problem P2-1 and a concave lower
v (t) = (Ω
bound C̄ji (vv ; v (t)) is derived as (37) in Appendix D for the MS i’s achievable rate constraint (20b). As
w opt
Ωopt
for the sum-rate maximization, we extract the feasible beamforming vectors {w
ji } from {Ω
ji } resulting
from the solution v̂v (vv (t)) with the standard rank-reduction method coupled with scaling method [38]. Its
convergence is established by the property of SCA method and monotonicity of the rate and the EFIM
functions with respect to the scaling factor s as discussed in Section IV-B.

19

V. K RONECKER B EAMFORMING D ESIGN
In this section, we investigate a reduced-complexity beamforming design based on the parameterization
of the Kronecker channel model (4). We tackle both the sum-rate maximization (P1) and the sumSPE minimization (P2). To this end, we apply a similar approach to that developed for unstructured
beamforming with the caveat that we tackle alternately the optimization of the azimuth and elevation
components by leveraging the decomposition of the FD channel model. The solution will also require to
estimate the azimuth and elevation components of the channel, which will be done by using the solution
called “Nearest Kronecker product” [39], [40].
Based on the decomposition of beamforming vectors w ji in (9), we start by writing the beamforming
covariance matrices Ω ji in the form Ω ji = Ω A,ji ⊗ Ω E,ji with the rank-1 azimuth covariance matrices
H
Ω A,ji = w A,jiw H
A,ji and the elevation covariance matrices Ω E,ji = w E,jiw E,ji for all j ∈ NB and i ∈ NM .

As mentioned, the azimuth components ĝg A,ji and elevation components ĝg E,ji of the estimated channel
ĝg ji = ĝg A,ji ⊗ ĝg E,ji are extracted by using the “Nearest Kronecker product” scheme [39], [40]. This
scheme finds the solution of the minimization problem minĝg A,ji ,ĝg E,ji kĝg ji − ĝg A,ji ⊗ ĝg E,ji kF based on the
singular value decomposition. Accordingly, the azimuth component and the elevation component are
computed as ĝg A,ji =

√

σji,1u ji,1 and ĝg E,ji =

√

σji,1v ji,1 , respectively, where when defining the NA,j × NE,j

Gji = U jiΣ jiV ji whose elements are taken columnwise from ĝg ji ; σji,1 is the largest singular
matrix Ĝ
value component of the matrix Σ ji ; and u ji,1 and v ji,1 are corresponding left and right singular vectors
of the matrix U ji and V ji , respectively. With the estimated channel components ĝg A,ji and ĝg E,ji , the
(k)
(k)
(k)
(k)
Ωjk ) in (11) is decomposed into ξji
Ωjk ) = (1 − Dj )2 ξA,ji
ΩA,jk )ξE,ji
ΩE,jk ), where
effective power ξji (Ω
(Ω
(Ω
(Ω
(k)
(k)
ΩA,jk ) = ĝg H
g A,ji and ξE,ji
ΩE,jk ) = ĝg H
g E,ji .
ξA,ji (Ω
(Ω
A,jiΩ A,jkĝ
E,jiΩ E,jkĝ

Ωji (t+1)}.
At each outer tth iteration, two inner loops are employed in order to obtain the next iterate {Ω
ΩA,ji (t + 1)} for fixed elevation covariance
The first is used to optimize the azimuth covariance matrices {Ω
ΩE,ji (t)}; while the second is used for optimizing the elevation covariance matrices {Ω
ΩE,ji (t +
matrices {Ω
ΩA,ji (t + 1)}. For each inner loop, the SCA-based approach
1)} with fixed azimuth covariance matrices {Ω

20

(a)

(b)

Fig. 3. (a) Set-up and (b) simulation parameters for the numerical results.

detailed in Algorithm 1 is applied to optimize the azimuth or elevation covariance matrices separately
κ, χ ).
along with the auxiliary variables (κ
VI. N UMERICAL R ESULTS
In this section, we evaluate the performance of the unstructured and Kronecker beamforming design
proposed in Section IV and Section V, respectively, in terms of the average sum-rate and the square
root of average sum-SPE via Monte Carlo simulations. As illustrated in Fig. 3(a), we consider a network
in which NB = 4 BSs are placed at the vertices of square region of side length ∆x = ∆y m at height
zB,j = 50 m for all BSs j ∈ NB , while the MSs are randomly uniformly distributed within a ∆x ×∆y ×∆z
cube. We consider both the maximization of the sum-rate under SPE constraints in problem P1 and the
minimization of the sum-SPE under rate constraints in problem P2. Furthermore, for reference, we include
the upper bound of the solution of problem P1 obtained by removing the localization accuracy constraints
(19b), as well as the lower bound on the optimal solution of problem P2 obtained by removing the data
rate constraints (20b). For both bounds, we consider the performance under unstructured beamforming.
Unless stated otherwise, we consider the Kronecker channel model as described in Section II-A2.
The parameters are summarized in Fig. 3(b) and described in this paragraph. We assume that each BS
j ∈ NB is equipped with Nj = 8 antennas consisting of NA,j = 2 horizontal antennas and NE,j = 4

21

(a)

(b)

Fig. 4. Average sum-rate for P1 and square root of average sum-SPE for P2 as a function of number B of output levels of ADC/DACs
at BSs with (a) correlated channel model and (b) Kronecker channel model (NM = 2 MSs randomly and uniformly distributed within the
cube area).

vertical antennas. We also set ∆x = ∆y = 500 m and ∆z = 50 m, and the reference distance d0 in
(2) is 0.01 so that the path loss at a distance of 100 m is ζ = −60 dB with path loss exponent η = 3.
Moreover, we assume a noise level of N0 = −174 dBm/Hz, the propagation speed c = 3 × 108 m/s, an
effective bandwidth of β = 100 MHz, and training and data phases with n̄p = np = 100 and nd = 200
symbols. The wavelength is chosen to be λ = 0.1667 m, which corresponds to carrier frequency 1.8 GHz,
and the azimuth and elevation angular perturbation are set to σφ = σθ = π/36. The total uplink and
downlink transmit powers are assumed to be NM P̄ = P = 30 dBm and the number of output levels at the
ADC/DACs at all the BSs are equal to B = Bj = 8 for j ∈ NB . Identical requirements for localization
accuracy ρmax = ρmax
= (0.01∆x )2 in problem P1 and data rate C min = Cimin = 0.15 in problem P2 are
i
applied to all MSs.
In Fig. 4, we investigate the effect of the resolution B of the ADC/DACs at the BSs with NM = 2
MSs in correlated channel model (Fig. 4(a)) and Kronecker channel model (Fig. 4(b)) by showing the
rate and localization performance obtained from problems (19) and (20), respectively. In general, as
the resolution B increases, both data rate and localization accuracy are improved. More interestingly, the

22

Fig. 5. Average sum-rate for P1 and square root of average sum-SPE for P2 as a function of number NE of vertical antennas at BSs
(Kronecker channel model; NM = 2 MSs randomly and uniformly distributed within the cube area).

proposed schemes with unstructured beamforming are able to accommodate localization or rate constraints
with only minor performance degradations as compared to the reference case without such constraints.
In contrast, reduced-complexity Kronecker solutions present a performance loss, which tends to decrease
with B. This is also observed to be slightly increased in correlated channel model as compared to the
Kronecker channel model. For instance, with B = 6 bits, the performance loss of Kronecker beamforming
over unstructured beamforming is 4.9% under the correlated channel model and 4.7% under the Kronecker
channel model in terms of data rate optimization in P1, while it amounts to 31.7% under the correlated
channel model and 24.9% loss under the Kronecker channel model in terms of localization accuracy in
P2.
Next, we consider the performance as a function of the number NE of vertical antennas in Fig. 5 with
the same setting with Fig. 4 but with B = 8 bits. Increasing the number NE of vertical antennas is seen to
provide enhanced data rate and localization performance due to the larger available number of degrees of
freedom. The figure confirms the main conclusions discussed above with regard to the effectiveness of the
proposed optimal strategy. Furthermore, a larger NE is also observed to reduce the performance loss of

23

Fig. 6. Average sum-rate for P1 and square root of average sum-SPE for P2 as a function of number NM of MSs (Kronecker channel
model; NM MSs randomly and uniformly distributed within the cube area).

Kronecker beamforming as compared to unstructured beamforming for rate optimization under localization
constraints. This reduction is instead less noticeable for localization optimization under rate constraints.
This can be explained by the fact that ensuring data rate requirements calls for the management of interMS interference, which is instead not an issue for positioning. In fact, localization accuracy is improved as
long as the overall received power on all beams is increased. This suggests that guaranteeing minimal rate
constraints requires a larger number of degrees of freedom than ensuring localization accuracy constraints.
Fig. 6 investigates the same performance criteria as a function of the number NM of MSs. First, it
is observed that a larger NM allows the achievement of an enhanced data rate performance, since the
presence of more user offers an increased multiuser diversity that allows the transmission to MSs that
have favorable channel conditions. In contrast, increasing NM yields a degraded performance in terms of
localization accuracy. This is because each MS tends to receive less power for a larger NM due to the
need to beamform to individual users in order to satisfy the rate constraints.

24

VII. C ONCLUDING R EMARKS
In this paper, we have investigated the problem of beamforming design for a TDD-based downlink
system with FD-MIMO BSs equipped with limited-resolution front ends. Unlike existing works, we
considered simultaneously the criteria of sum-rate data transmission and sum-SPE localization accuracy,
under the assumption that localization is performed at the MSs based on downlink signals. Two types
of beamforming designs are studied, namely unstructured beamforming and a low-complexity Kronecker
beamforming solution, where the latter decomposes the beamforming vector for each BS into separate
azimuth and elevation components. Extensive numerical results suggest that the proposed schemes can
effectively ensure both data and localization performance criteria with a limited performance loss as
compared to the case in which either data communication or positioning is considered. The interplay
between rate and localization accuracy is seen to arise from their different requirements on beamforming
design: While inter-MS interference management is required for the enhanced data rate in beamforming
design, positioning only calls for the maximization of the power received at the users, irrespective of
whether it is useful signal or interference. Optimized low-complexity Kronecker beamforming solutions
are observed to achieve limited performance loss when the front-end resolution is large enough, the number
of transmit antennas is large and the number of users is small. Among open issues left for future work,
we mention the performance comparison of downlink positioning and uplink positioning in asynchronous
scenarios.

A PPENDIX A
D ERIVATIONS OF γ1 , . . . , γ7

IN

(3)

In this appendix, by following [12], the definition of the γ1 , . . . , γ7 is recalled for the channel model
(3). The small-scale fading channel vector h ji can be written as h ji = ejςji a ji ⊗ b ji , where the phase
shift ςji of transmit path between BS j and MS i is assumed to be uniformly distributed in [0, 2π],
and the steering vectors for y-axis and z-axis are given as a ji = [1 e−j

2πdA
λ

cos(φji +∆φji ) sin(θji +∆θji )

···

25

e−j(NA,j −1)

2πdA
λ

] and b ji = [1 e−j

cos(φji +∆φji ) sin(θji +∆θji ) T

2πdE
λ

cos(θji +∆θji )

· · · e−j(NE,j −1)

2πdE
λ

cos(θji +∆θji ) T

] ,

respectively, for j ∈ NB and i ∈ NM , with the azimuth angular perturbation ∆φji ∼ N (0, σφ2 ) and elevation
angular perturbation ∆θji ∼ N (0, σθ2 ). As derived in [12], this yields the covariance matrix R hji in (3),
where γ1 = e
γ4 =

1
2

σθ 2π
λ

j

2πdE
λ

2



2πdE 2
1
(p−k)2 sin2 θji
(p−k) cos θji − 2 σθ λ

e

; γ2 =

2πdA
(q
λ

A
− l) sin θji ; γ3 = σθ 2πd
(q − l) cos θji ;
λ

dE dA (p − k)(q − l) sin(2θji ); γ5 = γ32 σφ2 sin2 φji + 1; γ6 = γ42 σφ2 sin2 φji + cos φji ; γ7 =

γ32 cos2 φji − γ42 σφ2 sin2 φji − 2γ4 cos φji .

A PPENDIX B
D ERIVATIONS OF (16)
Here, we derive the performance metric (16) for the localization accuracy of MS i. The MCRB is
the trace of the inverse of the average EFIM, where the average, in this work, is with respect to the
pilot signals [31], [32]. To evaluate the average EFIM, we start by defining the unknown parameter
W 1 ) · · · α TNB i
vector for MS i as ψ i = [ppTM,i α T1i σz̃21i (W

(1)
w j1 )}
W NB )]T , where α ji = [<{αji
(w
σz̃2N i (W
B

(NM )
(NM )
(1)
w jNM )}]T . Computing the average EFIM requires to evalw jNM )} ={αji
w j1 )} · · · <{αji
(w
(w
={αji (w

ψ i )2 ], which is averaged over the training sequences s p =
uate the quantity Ey pi ,ssp [(∂ ln f (yy pi |ssp , ψ i )/∂ψ
{sspji }j∈NB ,i∈NM , with s pji being the vector representation of spji (t) =

Pnp

l=1

spji (l)p(t − lTs ). This is made

p
difficult by the fact that the effective noise z̃ji
(l) in (12) is not independent of the useful signal and not

Gaussian.
To address this problem, we use the lower bound derived in [33], [34] that only requires the knowledge
of the first and second moment of the system output. Accordingly, we have the inequality
"
"
"
2 #
2 ##
∂ ln f (yy pi |ssp , ψ i )
∂ ln f (yy pi |ssp , ψ i )
Ey pi ,ssp
= Es p Ey pi |ssp
(31a)
ψi
ψi
∂ψ
∂ψ
"


#
µyip (ψ
ψ i) T
µyip (ψ
ψ i)
∂µ
∂µ
−1
ψ i)
W ), (31b)
 Es p
R yip (ψ
, J ψ i (W
ψi
ψi
∂ψ
∂ψ
where the inequality relationship in (31b) is obtained based on the Cauchy-Schwarz inequality; and
µyip (ψ
ψ i ) and R yip (ψ
ψ ) are the first and second output moments given by µyip (ψ
ψ i) =

R

y pi f (yy pi |ssp , ψ i ) dyy pi

26

ψ i) =
and R yip (ψ


T
R p
ψ i ) y pi − µ yip (ψ
ψ i ) f (yy pi |ssp , ψ i ) dyy pi , respectively. Note that we have emy i − µ yip (ψ

W j }j∈NB . To evaluate (31),
phasized the dependence of (31b) on the beamforming vectors W = {W
T

T

ψ NB i ]T with
ψ i , where ψ̃
ψ i = [ψ̃
ψ 1i · · · ψ̃
we relate the parameter vector ψ i to the larger parameter vector ψ̃
ψ ji = [τji α Tji σz̃2ji (W
W j )]T . When the MS is localizable, this mapping is a bijection [16], [29], [30], and
ψ̃
W ) in (31) as J ψ i (W
W ) = T J ψ̃ψ i (W
W )T
T T , with
this reparameterization allows us to rewrite the matrix J ψ i (W
ψ i and J ψ̃ψ i (W
W ) being defined in the same
T being the Jacobian matrix for the transformation from ψ i to ψ̃
W ) in (31b) with ψ̃
ψ i in lieu of ψ i . By applying the standard Schur complement condition
way of J ψ i (W
W ), we obtain the following lower bound on the submatrix of J ψ i (W
W ) in (17)
[41] to the matrix J ψ i (W
corresponding to the MS i’s position p M,i . The derivation of (17) from (31) follows using the same steps
as in [16], [29], [30], which yields the performance metric for the SPE in (15) of MS i as the MCRB in
(16).

A PPENDIX C
D ERIVATION OF P ROBLEM P1-2
In this appendix, we derive the SCA algorithm in Algorithm 1 for problem P1-1. We first observe that
the objective function in (23a) is the sum of rates Cji (vv ), each having the DC form
Cji (vv ) =


nd
fji+ (vv ) − fji− (vv ) ,
NB n

(32)

where




Nj

X

fji+ (vv ) = log2 

(k)
Ωjk ) + N0 +
ξji (Ω

k∈NM

X

R ∆gji Ω jk } +
(1 − Dj )2 tr{R

k∈NM

X


κjn tr E nR gji E Tn 

(33a)

n=1

and

fji− (vv ) = log2 

X

(k)
Ωjk ) + N0 +
ξji (Ω

k∈NM ,k6=,i

X
k∈NM

R ∆gji Ω jk } +
(1 − Dj )2 tr{R

Nj
X



κjn tr E nR gji E Tn  .(33b)

n=1

Then, given v (t) for the tth feasible iterate, we obtain the strongly concave surrogate function Ĉ(vv ; v (t))
in (26) by applying Lemma 1.

27

For the non-convex constraint (23d), we first define the function


Nj
X
X

R ∆gji Ω jk } +
gji (v1 , v 2 ) , χji 
κjn tr E nR gji E Tn 
(1 − Dj )2 tr{R
=

(34a)

n=1

k∈NM


1 2
1
(h1,ji (v1 ) + h2,ji (vv 2 ))2 −
h1,ji (v1 ) + h22,ji (vv 2 ) ,
2
2

for j ∈ NB and i ∈ NM , where h1,ji (v1 ) = χji ; h2,ji (vv 2 ) =

P

k∈NM (1

(34b)

R ∆gji Ω jk } +
− Dj )2 tr{R

PNj

n=1

κjn


Ωj1 , . . . , Ω jNM , κj1 , . . . , κjNj }. Given a feasible solution v (t), using
tr E nR gji E Tn ; v1 = χji and v 2 = {Ω
Lemma 2, we then derive the upper bound ḡji (v1 , v 2 ; v1 (t), v 2 (t)) ≥ gji (v1 , v 2 ), where the left-hand side is
defined as (29). Consequently, since the constraint (23d) can be written as gji (v1 , v 2 ) ≤
χji N0 , it can be convexified as ḡji (v1 , v 2 ; v1 (t), v 2 (t)) ≤

P

k∈NM

P

k∈NM

(k)
Ωjk )−
ξji (Ω

(k)
Ωjk )−χji N0 in (28b) while satisfying
ξji (Ω

the conditions of the SCA method (Lemma 2).

A PPENDIX D
D ERIVATION OF P ROBLEM P2-2
In this appendix, we provide technical details on the derivations of the SCA method in Algorithm 1
for the sum-SPE minimization problem. We start by relaxing the rank constraint on Ω and obtain the
equivalent problem:

!−1 
N
B


2
2
X 1
8π np β X
J
(P2-1): minimize
tr
(35a)
φ,θ (φji,m , θji,m )χji
Ω,κ
κ,χ
χ

Ns m=1 
c2
j=1
i∈NM


Nj
X
X
X
X

nd
(k)
Ωjk ) + N0 +
R ∆gji Ω jk } +
s.t.
log2 
ξji (Ω
(1 − Dj )2 tr{R
κjn tr E nR gji E Tn 
NB n j∈N
n=1
k∈NM
k∈NM
B


Nj
X
X
X

(k)
Ωjk ) + N0 +
R ∆gji Ω jk } +
− log2 
ξji (Ω
(1 − Dj )2 tr{R
κjn tr E nR gji E Tn  ≥ Cimin , i ∈ NM
Ns
X

k∈NM ,k6=,i

k∈NM

n=1

(35b)
(23c) - (23g).

(35c)

The problem P2-1 is not convex due to the non-convex constraints (35b) and (23d). In order to apply the
SCA method, we first recall the following lemma.

28

x) ≥ 0, where gi (x
x) has the DC
Lemma 3 ([43, Example 3]): Consider a non-convex constraint gi (x
x) = gi+ (x
x) − gi− (x
x), with gi+ (x
x) and gi− (x
x) being concave and continuously differenstructure, namely gi (x
x), a concave lower approximation ḡi (x
x; y ) ≤ gi (x
x) that guarantees
tiable. For any y in the domain of gi (x
the requirements [36, Assumption 3] of the SCA algorithm is given as
x; y ) , gi+ (x
x) − gi− (yy ) − ∇x gi− (yy )T (x
x − y ).
ḡi (x

(36)

Ω, κ , χ ) and v (t) = (Ω
Ω(t), κ (t), χ (t))
We define the set of optimization variables for problem P2-1 as v = (Ω
for the tth iterate within the feasible set of the problem P2-1. In addition, the non-convex constraints
(35b) can be written as

P

j∈NB

Cji (vv ) ≥ Cimin in (35b), where Cji (vv ) is a DC function with fji+ (vv ) and

+
−
fji− (vv ) in (32) as gji
(vv ) and gji
(vv ), respectively. As a result, using Lemma 3, a concave lower bound

C̄ji (vv ; v (t)) ≤ Cji (vv ) can be derived for use in the SCA algorithm as
C̄ji (vv ; v (t)) ,


nd
+
−
−
gji
(vv ) − gji
(vv (t)) − ∇v gji
(vv (t))T (vv − v (t)) ,
NB n

(37)

−
where ∇v gji
(vv (t))T (vv − v (t)) is calculated as (27). For the non-convex constraints (23d), a convex

upper bound ḡji (v1 , v 2 ; v1 (t), v 2 (t)) ≥ gji (v1 , v 2 ) in (29) is obtained using Lemma 2, yielding the convex
constraint (28b). Consequently, for a feasible v (t), we obtain problem P2-2 in (30).

R EFERENCES
[1] R. D. Taranto, S. Muppirisetty, R. Raulefs, D. Slock, T. Svensson, and H. Wymeersch, “Location-aware communications for 5G
networks: How location information can improve scalability, latency, and robustness of 5G,” IEEE Sig. Proc. Mag., vol. 31, no. 6, pp.
102–112, Nov. 2014.
[2] X. Lin, J. Bergman, F. Gunnarsson, O. Liberg, S. M. Razavi, H. S. Razaghi, H. Rydén, and Y. Sui, “Positioning for the internet of
things: A 3GPP perspective,” arXiv preprint arXiv:1705.04269v1, Apr. 2017.
[3] 5G-PPP,

“5G

empowering

vertical

industries,”

Feb.

2015,

[Online].

Available:

https://5g-ppp.eu/wp-

content/uploads/2016/02/BROCHURE 5PPP BAT2 PL.pdf.
[4] C. Hoymann, D. Astely, M. Stattin, G. Wikstrom, J. F. Cheng, A. Hoglund, M. Frenne, R. Blasco, J. Huschke, and F. Gunnarsson,
“LTE release 14 outlook,” IEEE Comm. Mag., vol. 54, no. 6, pp. 44–49, Jun. 2016.
[5] G. T. 36.897, “Study on elevation beamforming/full-dimension (FD) multiple input multiple output (MIMO) for LTE (Release 13).”

29

[6] Y.-H. Nam, B. L. Ng, K. Sayana, Y. Li, J. Zhang, Y. Kim, and J. Lee, “Full-dimension MIMO (FD-MIMO) for next generation cellular
technology,” IEEE Comm. Mag., vol. 51, no. 4, pp. 172–179, Jun. 2013.
[7] Y. Li, C. Tao, G. Seco-Granados, A. Mezghani, A. L. Swindlehurst, and L. Liu, “Channel estimation and performance analysis of
one-bit massive MIMO systems,” IEEE Trans. Sig. Proc., vol. 65, no. 15, pp. 4075–4089, Aug. 2017.
[8] S. Jacobsson, G. Durisi, M. Coldrey, T. Goldstein, and C. Studer, “Quantized precoding for massive MU-MIMO,” IEEE Trans. Comm.,
vol. PP, no. 99, Jul. 2017.
[9] B. Mondal, T. A. Thomas, E. Visotsky, F. W. Vook, A. Ghosh, Y.-H. Nam, Y. Li, J. Zhang, M. Zhang, Q. Luo, Y. Kakishima, and
K. Kitao, “3d channel model in 3GPP,” IEEE Comm. Mag., vol. 53, no. 3, pp. 16–23, Mar. 2015.
[10] A. Kammoun, H. Khanfir, Z. Altman, M. Debbah, and M. Kamoun, “Preliminary results on 3D channel modeling: From theory to
standardization,” IEEE Jour. Select. Areas in Comm., vol. 32, no. 6, pp. 1219–1229, June 2014.
[11] J. Kang, O. Simeone, J. Kang, and S. Shamiai, “Layered downlink precoding for C-RAN systems with full dimensional MIMO,” IEEE
Trans. Veh. Technol., vol. 66, no. 3, pp. 2170–2182, Mar. 2017.
[12] D. Ying, F. W. Vook, T. A. Thomas, D. J. Love, and A. Ghosh, “Kronecker product correlation model and limited feedback codebook
design in a 3D channel model,” in Proc. IEEE Int. Conf. on Comm., Jun. 2014, Sydney, Australia.
[13] C. G. Tsinos, S. Maleki, S. Chatzinotas, and B. Ottersten, “Hybrid analog-digital transceiver designs for cognitive radio millimiter
wave systems,” in Proc. of Asilomar Conf. on Sign., Syst. and Computers, Nov. 2016, pp. 1785–1789, Pacific Grove, CA, USA.
[14] S. Wang, K. D. Lee, S. G. Kim, and J. S. Kwak, “Enhance downlink positioning in WiMAX/16m.”
[15] J. J. Bussgang, “Crosscorrelation functions of amplitude-distorted Gaussian signals,” Res. Lab. Elec. Tech. Rep. 216, Mar. 1952,
cambridge, MA, USA.
[16] S. Jeong, O. Simeone, A. Haimovich, and J. Kang, “Beamforming design for joint localization and data transmission in distributed
antenna system,” IEEE Trans. Veh. Technol., vol. 64, no. 1, pp. 62–76, Jan. 2015.
[17] T. L. Marzetta, “How much training is required for multiuser MIMO?” in Proc. of Asilomar Conf. on Sign., Syst. and Computers, Oct.
2006, Pacific Grove, CA, USA.
[18] S. Shahsavari, P. Hassanzadeh, A. Ashikhmin, and E. Erkip, “Sectoring in multi-cell massive MIMO systems,” arXiv preprint
arXiv:1707.09070, Jul. 2017.
[19] A. Mezghani and J. A. Nossek, “Capacity lower bound of MIMO channels with output quantization and correlated noise,” in Proc.
IEEE Int. Symp. Info. Th., Jul. 2012, Cambridge, MA, USA.
[20] M. Stein, F. Wendler, A. Mezghani, and J. Nossek, “Quantization-loss reduction for signal parameter estimation,” in Proc. IEEE Int.
Conf. Acoust., Speech and Sig. Proc., May 2013, Vancouver, Canada.
[21] A. Kakkavas, J. Munir, A. Mezghani, H. Brunner, and J. Nossek, “Weighted sum rate maximization for multi-user MISO systems with
low resolution digital to analog converters,” in Proc. International ITG Workshop on Smart Antennas, Mar. 2016, Munich, Germany.
[22] A. F. Molisch, Wireless communications.

Wiley, 2010.

[23] A. Gershman and N. Sidiropoulos, Space-time processing for MIMO communications.

New york: John Wiley & Sons, 2005.

30

[24] J. Max, “Quantizing for minimum distortion,” IRE Transactions on Information Theory, vol. 6, no. 1, pp. 7–12, Mar. 1960.
[25] B. Hassibi and B. M. Hochwald, “How much training is needed in multiple-antenna wireless links?” IEEE Trans. Info. Th., vol. 49,
no. 4, Apr. 2003.
[26] A. E. Gamal and Y.-H. Kim, Network Information Theory.

Cambridge University Press, 2011.

[27] M. Medard, “The effect upon channel capacity in wireless communications of perfect and imperfect knowledge of the channel,” IEEE
Trans. Info. Th., vol. 46, no. 3, pp. 933–946, May 2000.
[28] S. M. Kay, Fundamentals of Signal Processing-Estimation Theory.

Englandwood Cliffs, NJ: Prentice Hall, 1993.

[29] S. Jeong, O. Simeone, A. Haimovich, and J. Kang, “Optimal fronthaul quantization for cloud radio positioning,” IEEE Trans. Veh.
Technol., vol. 65, no. 4, pp. 2763–2768, Apr. 2016.
[30] ——, “Positioning via direct localization in C-RAN systems,” IET communications, vol. 10, no. 16, pp. 2238–2244, Nov. 2016.
[31] A. N. D’Andrea, U. Mengali, and R. Reggiannini, “The modified Cramèr-Rao bound and its application to synchronization problems,”
IEEE Trans. Comm., vol. 42, no. 234, pp. 1391–1399, Feb./Mar./Apr. 1994.
[32] F. Gini, R. Reggiannini, and U. Mengali, “The modified Cramèr-Rao bound in vector parameter estimation,” IEEE Trans. Comm.,
vol. 42, no. 234, pp. 1391–1399, Feb./Mar./Apr. 1994.
[33] M. Stein, A. Mezghani, and J. A. Nossek, “A lower bound for the fisher information measure,” IEEE Sig. Proc. Lett., vol. 21, no. 7,
pp. 796–799, Jul. 2014.
[34] M. Stein and J. A. Nossek, “A pessimistic approximation for the fisher information measure,” IEEE Trans. Sig. Proc., vol. 65, no. 2,
pp. 386–396, Jan. 2017.
[35] A. Shapiro, D. Dentcheva, and A. Ruszczyski, Lectures on stochastic programming: modeling and theory. Society for Industrial and
Applied Mathematics (SIAM), 2009.
[36] G. Scutari, F. Facchinei, L. Lampariello, and P. Song, “Parallel and distributed methods for nonconvex optimizationpart I: Theory,”
arXiv preprint arXiv:1410.4754v2, Jan. 2016.
[37] G. Scutari, F. Facchinei, L. Lampariello, P. Song, and S. Sardellitti, “Parallel and distributed methods for nonconvex optimizationpart
II: Applications,” arXiv preprint arXiv:1601.04059v1, Jan. 2016.
[38] Z.-Q. Luo, W.-K. Ma, A.-C. So, Y. Ye, and S. Zhang, “Semidefinite relaxation of quadratic optimization problems,” IEEE Signal Proc.
Magazine, vol. 27, no. 3, pp. 20–34, May 2010.
[39] G. H. Golub and C. V. Loan, Matrix computations.

The John Hopkins University Press, 2012.

[40] C. V. Loan and N. Pitsianis, Approximation with Kronecker products.
[41] F. Zhang, The Schur complement and its applications.

Cornell University, 1992.

Springer, 2005, vol. 4.

