# Pilot Reuse Among D2D Users in D2D Underlaid Massive MIMO Systems

Hao Xu, _Student Member, IEEE_ 0, Wei Xu, _Senior Member, IEEE_ 0, Zhaohui
Yang, _Student Member, IEEE_ 0, Jianfeng Shi, _Student Member, IEEE_ 0, and
Ming Chen, _Member, IEEE_ 0 Copyright (c) 2015 IEEE. Personal use of this
material is permitted. However, permission to use this material for any other
purposes must be obtained from the IEEE by sending a request to pubs-
<EMAIL_ADDRESS>This work was in part supported by the NSFC (Nos.
61372106, 61471114, & 61221002), NSTMP under 2016ZX03001016-003, the Six
Talent Peaks project in Jiangsu Province under GDZB-005, Science and
Technology Project of Guangdong Province under Grant 2014B010119001, the
Scholarship from the China Scholarship Council (No. 201606090039), Program
Sponsored for Scientific Innovation Research of College Graduate in Jiangsu
Province under Grant KYLX16_0221, and the Scientific Research Foundation of
Graduate School of Southeast University under Grant YBJJ1651. _(Corresponding
author: Hao Xu, Wei Xu.)_ H. Xu, W. Xu, Z. Yang, J. Shi and M. Chen are with
the National Mobile Communications Research Laboratory, Southeast University,
Nanjing 210096, China (Email: {xuhao2013, wxu, yangzhaohui, shijianfeng and
chenming}@seu.edu.cn).

###### Abstract

In a device-to-device (D2D) underlaid massive MIMO system, D2D transmitters
reuse the uplink spectrum of cellular users (CUs), leading to cochannel
interference. To decrease pilot overhead, we assume pilot reuse (PR) among D2D
pairs. We first derive the minimum-mean-square-error (MMSE) estimation of all
channels and give a lower bound on the ergodic achievable rate of both
cellular and D2D links. To mitigate pilot contamination caused by PR, we then
propose a pilot scheduling and pilot power control algorithm based on the
criterion of minimizing the sum mean-square-error (MSE) of channel estimation
of D2D links. We show that, with an appropriate PR ratio and a well designed
pilot scheduling scheme, each D2D transmitter could transmit its pilot with
maximum power. In addition, we also maximize the sum rate of all D2D links
while guaranteeing the quality of service (QoS) of CUs, and develop an
iterative algorithm to obtain a suboptimal solution. Simulation results show
that the effect of pilot contamination can be greatly decreased by the
proposed pilot scheduling algorithm, and the PR scheme provides significant
performance gains over the conventional orthogonal training scheme in terms of
system spectral efficiency.

##  I Introduction

With the increasing demand on broadband wireless communications, the problem
of spectrum insufficiency has become a major factor limiting the wireless
system performance [1]. Massive multiple-input multiple-output (MIMO)
transmission was proposed in [2] and has triggered considerable research
interest recently due to its great gains in spectral efficiency (SE) and
energy efficiency (EE) [3, 4, 5]. Besides, device-to-device (D2D)
communication has also been proven promising in enhancing the SE of the
traditional cellular systems and has drawn great attention recently [6, 7, 8].
Different from the conventional cellular communication where all traffic is
routed via base station (BS), D2D communication allows two closely located
users to communicate directly, and thus have distinct advantages such as high
SE, short packet delay, low energy consumption and increased safety.

There has been extensive research on design and analysis of massive MIMO
systems [9, 10, 11, 12, 13]. In [9], the uplink capacity bounds were derived
under both perfect and imperfect channel state information (CSI), and the
tradeoff between SE and EE was studied. Ref. [10] compared two most prominent
linear precoders with respect to (w.r.t.) SE and radiated EE in a massive MIMO
system. Unlike [9] and [10], which considered simplified single-cell
scenarios, [11, 12, 13] studied multi-cell massive MIMO systems. As for
underlaid D2D communication, a great challenge to the existing cellular
architecture is the cochannel interference due to spectrum reuse. There has
been a lot of literature working on interference mitigation for D2D underlaid
systems [14, 15, 16, 17, 18, 19, 20]. In [14, 15, 16, 17], resource allocation
and power control algorithms were proposed to maximize the SE of D2D users
(DUs), and in [18, 19, 20], extended algorithms were carried out to maximize
the EE of DUs.

Though massive MIMO and D2D communication have been widely studied, only a few
papers investigated the interplay between massive MIMO and D2D communication
[21, 22]. In [22], the SE of cellular and D2D links was investigated under
both perfect and imperfect CSI, but the overhead for acquiring CSI was not
considered. In massive MIMO systems, orthogonal pilots are transmitted by
cellular users (CUs) to obtain CSI. When D2D communication is introduced and
orthogonal pilots are used at each D2D transmitter (D2D-Tx) for channel
estimation, the pilot overhead is large which will significantly affect the
system performance. Furthermore, as a multi-user transmission strategy,
massive MIMO is designed to support multiple users transmitting on the same
time-frequency block. Though D2D-to-cellular interference can be greatly
reduced by a large antenna array at BS, cellular-to-D2D interference still
persists and may be worse than a conventional single-input single-output
(SISO) D2D underlaid system.

In order to shorten pilot overhead, an effective strategy is allowing
orthogonal pilots to be reused among different users. Most of the existing
works with pilot reuse (PR) mainly focus on multiple-cell scenarios, i.e.,
mobile users in the same cell use orthogonal pilots, and users in different
cells reuse the same set of pilots [23, 24, 25, 26, 27]. To the best of the
authors’ knowledge, only a few works have considered the strategy of PR within
a cell [28, 29, 30, 31]. In [28], the authors analyzed the feasibility of PR
over spatially correlated massive MIMO channels with constrained channel
angular spreads. Authors of [29] and [30] allowed D2D-Txs to reuse the pilots
of CUs and proposed an interference-aided minimum-mean-square-error (MMSE)
detector to suppress the D2D-to-cellular interference. [31] also studied a D2D
underlaid massive MIMO system with PR, but the performance of CUs was left out
of consideration for simplicity. In contrast to these existing works, our work
analyzes the achievable rate of both cellular and D2D links under PR, and
proposes pilot scheduling as well as power control algorithms to optimize the
system performance. The main contributions of this paper are summarized as
follows:

$\bullet$ We assume that CUs use orthogonal pilots while all D2D-Txs reuse
another set of pilots for channel estimation. The motivation of PR stems from
that D2D pairs usually locate dispersively and use low power for short-
distance transmission. Hence, letting several D2D pairs which are far away
from each other use the same pilot for channel estimation would cause
endurable pilot contamination. Under PR, we first derive the expression of
MMSE estimate of all channels. With the obtained imperfect CSI, all receivers
apply the partial zero forcing (PZF) receive filters studied in [32] for
signal detection. Then, we derive the effective signal-to-interference-plus-
noise ratio (SINR) and a lower bound on the ergodic achievable rate of each
user.

$\bullet$ Different from the estimation of cellular channel vectors which is
only affected by noise, the channel estimation of D2D links experiences effect
from both noise and pilot contamination due to PR. To mitigate pilot
contamination, we develop a pilot scheduling and pilot power control algorithm
under the criterion of minimizing the sum mean-square-error (MSE) of channel
estimation of D2D links. We first show that with an appropriate number of
orthogonal pilots available for DUs and a well designed pilot scheduling
scheme, each D2D-Tx should transmit its pilot using the maximum power. Then,
we develop a heuristic pilot scheduling scheme to allocate pilots to DUs, and
show that the sum MSE of channel estimation of D2D links can be decreased
significantly.

$\bullet$ We study the sum SE maximization for D2D links while guaranteeing
the quality of service (QoS) of CUs. Efficient power control algorithms often
play an important role in reducing cochannel interference and reaping the
potential benefits of D2D communication. However, these algorithms are usually
carried out based on the knowledge of instantaneous CSI of all links [14, 15,
16, 17]. Apart from the computational complexity, such algorithms require BS
to gather instantaneous CSI of all links, which is difficult for
implementation. Therefore, in this paper, we consider performing the power
control algorithm periodically at a coarser frame level granularity based on
large-scale fading coefficients which vary slowly. Simulation results show
that the proposed algorithm converges rapidly and obtains a much higher sum SE
of D2D links compared to the typical orthogonal training scheme.

Note that since PR among D2D pairs in a D2D underlaid massive MIMO system has
been less well researched, we consider a simplified single-input multiple-
output (SIMO) transmission for D2D communication as [22], and mainly focus on
analyzing the effect of PR on the system performance. As for MIMO
transmission, similar results can be obtained by simply modifying the analysis
and optimization of this manuscript. On the other hand, our analysis in this
paper focuses on a single-cell scenario for the sake of clarity. Regarding the
multi-cell massive MIMO system, there has been a lot of literature working on
mitigating pilot contamination [33, 34, 35]. As a result, for a multi-cell D2D
underlaid massive MIMO system, where PR among D2D pairs persists, we can first
use the algorithms developed in [33, 34, 35] to allocate pilots to CUs if CSI
can be exchanged among cells, and then straightforwardly extend the proposed
algorithms to improve system performance.

In this paper, we follow the common notations. $\mathbb{N}$, $\mathbb{R}$ and
$\mathbb{C}$ denote the set of natural numbers, the real space and the complex
space, respectively. The boldface upper (lower) case letters are used to
denote matrices (vectors). ${\bm{I}}_{M}$ stands for the $M\times M$
dimensional identity matrix and $\bm{0}$ denotes the all-zero vector or
matrix. “ $\setminus$ ” represents the set subtraction operation. Superscript
$(\cdot)^{H}$ denotes the conjugated-transpose operation and
${\mathbb{E}}\\{\cdot\\}$ denotes the expectation operation. We use
$\left\|{\bm{a}}\right\|_{2}$ to denote the Euclidean norm of $\bm{a}$.
${\bm{a}}\succeq\bm{0}({\bm{a}}\succ\bm{0})$ means that each element in
${\bm{a}}$ is positive (nonnegative).

The rest of this paper is organized as follows. In Section II, a D2D underlaid
massive MIMO system is introduced. In Section III, we present the MMSE
estimate of all channels under PR and show how PR affects the channel
estimation. The achievable rate of both cellular and D2D links is analyzed in
Section IV. In Section V, we aim to minimize the sum MSE of channel estimation
of D2D links and maximize the sum SE of all D2D links. Finally, numerical
verifications are presented in Section VI before concluding remarks in Section
VII.

##  II System Model

Consider the uplink of a D2D underlaid massive MIMO system with one BS, $N$
CUs and $K$ D2D pairs. The set of CUs and D2D pairs are denoted by ${\cal
N}=\\{1,\cdots,N\\}$111Here we misuse the notation ${\cal N}$ while avoiding
possible ambiguity with the ${\cal N}$ in the complex normal distribution sign
${\cal CN}$. and ${\cal{K}}=\\{1,\cdots,K\\}$, respectively. The BS is
equipped with $B$ antennas and each CU has one transmit antenna. As for the
D2D communication, we assume SIMO transmission, i.e., each D2D-Tx is equipped
with one antenna and each D2D receiver (D2D-Rx) is equipped with $M$ antennas
as in [22]. In this system, all transmitters use the same time-frequency
resource block to transmit signals, leading to cochannel interference. The
$B\times 1$ dimensional received data vector at BS is

${\bm{y}}^{(\text{c})}\\!=\\!\sum\limits_{n=1}^{N}\sqrt{q_{{\text{s}},n}u_{n}^{(\text{c})}}{\bm{h}}_{n}^{(\text{c})}x_{n}^{(\text{c})}\\!+\\!\sum\limits_{i=1}^{K}\sqrt{p_{{\text{s}},i}u_{i}^{(\text{d})}}{\bm{h}}_{i}^{(\text{d})}x_{i}^{(\text{d})}\\!+{\bm{z}},$
(1)

where $q_{{\text{s}},n}$ is the data transmit power of CU $n$, and
$x_{n}^{(\text{c})}$ is the zero-mean unit-variance data symbol of CU $n$.
$u_{n}^{(\text{c})}$ is the real-valued large-scale fading coefficient from CU
$n$ to BS and is assumed to be known as a priori.
${\bm{h}}_{n}^{(\text{c})}\sim{\cal CN}({\bm{0}},{\bm{I}}_{B})$ denotes the
fast fading vector channel from CU $n$ to BS. $p_{{\text{s}},i}$,
$x_{i}^{(\text{d})}$, $u_{i}^{(\text{d})}$ and ${\bm{h}}_{i}^{(\text{d})}$ are
similarly defined for D2D-Tx $i$. ${\bm{z}}\in{\mathbb{C}}^{B}$ is the complex
Gaussian noise at BS with covariance $N_{0}{\bm{I}}_{B}$.

Analogously, the $M\times 1$ dimensional received data vector at D2D-Rx $k$ is
given by

${\bm{y}}_{k}^{(\text{d})}\\!\\!=\\!\sum\limits_{i=1}^{K}\\!\sqrt{p_{{\text{s}},i}v_{ik}^{(\text{d})}}{\bm{g}}_{ik}^{(\text{d})}x_{i}^{(\text{d})}\\!+\\!\sum\limits_{n=1}^{N}\\!\sqrt{q_{{\text{s}},n}v_{nk}^{(\text{c})}}{\bm{g}}_{nk}^{(\text{c})}x_{n}^{(\text{c})}\\!+{\bm{n}}_{k},$
(2)

where $v_{ik}^{(\text{d})}$ and ${\bm{g}}_{ik}^{(\text{d})}\sim{\cal
CN}({\bm{0}},{\bm{I}}_{M})$ denote the real-valued large-scale fading
coefficient and the fast fading vector channel from D2D-Tx $i$ to D2D-Rx $k$,
respectively. $v_{nk}^{(\text{c})}$ and ${\bm{g}}_{nk}^{(\text{c})}$ are
similarly defined for the link from CU $n$ to D2D-Rx $k$.
${\bm{n}}_{k}\in{\mathbb{C}}^{M}$ is the complex Gaussian noise at D2D-Rx $k$
with covariance $N_{0}{\bm{I}}_{M}$.

##  III Channel Estimation

Orthogonal pilots are usually adopted to obtain the CSI of all links. In a D2D
underlaid system, to reduce pilot overhead, we assume that CUs use orthogonal
pilots while all D2D-Txs reuse another set of pilots for channel estimation.
Denote
${\bm{\Omega}}=\left({\bm{\omega}}_{1},\cdots,{\bm{\omega}}_{N},{\bm{\omega}}_{N+1},\cdots,{\bm{\omega}}_{\tau}\right)\in$
${\mathbb{C}}^{\tau\times\tau}$ as the pilot matrix with orthogonal column
vectors (i.e., ${\bm{\Omega}}^{H}{\bm{\Omega}}\\!=\\!{\bm{I}}_{\tau}$). $\tau$
($N<\tau\leq N+K$) is the length of the pilots and is also the number of
pilots available for channel estimation (this is the smallest amount of pilots
that are required). Then, without loss of generality, we assume that pilot
${\bm{\omega}}_{n}$ ($1\leq n\leq N$) is allocated to CU $n$ and the remaining
pilots $\\{{\bm{\omega}}_{N+1},\cdots,{\bm{\omega}}_{\tau}\\}$ are reused
among all D2D pairs.

###  III-A Channel Estimation at BS

Similar as the uplink data transmission in (1), the $B\times\tau$ dimensional
received signal matrix of pilot transmission at BS is

${\bm{Y}}^{(\text{c})}\\!=\\!\sum\limits_{n=1}^{N}\\!\sqrt{q_{{\text{p}},n}u_{n}^{(\text{c})}}{\bm{h}}_{n}^{(\text{c})}{\bm{\omega}}_{n}^{H}\\!+\\!\sum\limits_{i=1}^{K}\\!\sqrt{p_{{\text{p}},i}u_{i}^{(\text{d})}}{\bm{h}}_{i}^{(\text{d})}{\bm{\lambda}}_{i}^{H}\\!+\\!{\bm{Z}},$
(3)

where $q_{{\text{p}},n}$ and $p_{{\text{p}},i}$ denote the pilot transmit
powers of CU $n$ and D2D-Tx $i$, respectively.
${\bm{\lambda}}_{i}\in\\{{\bm{\omega}}_{N+1},\cdots,{\bm{\omega}}_{\tau}\\}$
is the pilot allocated to D2D pair $i$. ${\bm{Z}}$ is the noise matrix which
consists of independently and identically distributed (i.i.d.) Gaussian
elements with zero mean and variance $N_{0}$. Then, the MMSE estimate of
${\bm{h}}_{n}^{(\text{c})}$ is given by [36]

${\hat{\bm{h}}_{n}^{(\text{c})}}=\frac{\sqrt{q_{{\text{p}},n}u_{n}^{(\text{c})}}}{q_{{\text{p}},n}u_{n}^{(\text{c})}+N_{0}}{\bm{Y}}^{(\text{c})}{\bm{\omega}}_{n},\forall
n\in{\cal N}.$ (4)

Given the channel estimate vector ${\hat{\bm{h}}_{n}^{(\text{c})}}$, we can
express the true channel vector ${\bm{h}}_{n}^{(\text{c})}$ as
${\bm{h}}_{n}^{(\text{c})}={\hat{\bm{h}}_{n}^{(\text{c})}}+{\tilde{\bm{h}}_{n}^{(\text{c})}}$,
where the error vector ${\tilde{\bm{h}}_{n}^{(\text{c})}}$ represents the CSI
uncertainty. Due to the property of MMSE estimation [36],
${\hat{\bm{h}}_{n}^{(\text{c})}}$ is statistically independent of
${\tilde{\bm{h}}_{n}^{(\text{c})}}$, and they follow

${\hat{\bm{h}}_{n}^{(\text{c})}}\sim{\cal
CN}({\bm{0}},\delta_{n}^{(\text{c})}{\bm{I}}_{B}),{\tilde{\bm{h}}_{n}^{(\text{c})}}\sim{\cal
CN}({\bm{0}},\varepsilon_{n}^{(\text{c})}{\bm{I}}_{B}),$ (5)

where

$\delta_{n}^{(\text{c})}=\frac{q_{{\text{p}},n}u_{n}^{(\text{c})}}{q_{{\text{p}},n}u_{n}^{(\text{c})}+N_{0}},\varepsilon_{n}^{(\text{c})}=1-\delta_{n}^{(\text{c})}.$
(6)

We then analogously derive the MMSE estimate of ${\bm{h}}_{k}^{(\text{d})}$ as
follows

${\hat{\bm{h}}_{k}^{(\text{d})}}=\frac{\sqrt{p_{{\text{p}},k}u_{k}^{(\text{d})}}}{\sum\limits_{i\in{\cal
X}_{k}}p_{{\text{p}},i}u_{i}^{(\text{d})}+N_{0}}{\bm{Y}}^{(\text{c})}{\bm{\lambda}}_{k},\forall
k\in{\cal K},$ (7)

where ${\cal X}_{k}$ is the set of all D2D pairs using the same pilot as D2D
pair $k$. Denote
${\bm{h}}_{k}^{(\text{d})}={\hat{\bm{h}}_{k}^{(\text{d})}}+{\tilde{\bm{h}}_{k}^{(\text{d})}}$.
Then, ${\hat{\bm{h}}_{k}^{(\text{d})}}$ and
${\tilde{\bm{h}}_{k}^{(\text{d})}}$ are statistically independent satisfying

${\hat{\bm{h}}_{k}^{(\text{d})}}\sim{\cal
CN}({\bm{0}},\delta_{k}^{(\text{d})}{\bm{I}}_{B}),{\tilde{\bm{h}}_{k}^{(\text{d})}}\sim{\cal
CN}({\bm{0}},\varepsilon_{k}^{(\text{d})}{\bm{I}}_{B}),$ (8)

where

$\delta_{k}^{(\text{d})}=\frac{p_{{\text{p}},k}u_{k}^{(\text{d})}}{\sum\limits_{i\in{\cal
X}_{k}}p_{{\text{p}},i}u_{i}^{(\text{d})}+N_{0}},\varepsilon_{k}^{(\text{d})}=1-\delta_{k}^{(\text{d})}.$
(9)

###  III-B Channel Estimation at D2D-Rxs

The $M\times\tau$ dimensional received pilot signal matrix at D2D-Rx $k$ can
be written as

${\bm{Y}}_{k}^{(\text{d})}\\!=\\!\sum\limits_{n=1}^{N}\sqrt{q_{{\text{p}},n}v_{nk}^{(\text{c})}}{\bm{g}}_{nk}^{(\text{c})}{\bm{\omega}}_{n}^{H}\\!+\\!\sum\limits_{i=1}^{K}\sqrt{p_{{\text{p}},i}v_{ik}^{(\text{d})}}{\bm{g}}_{ik}^{(\text{d})}{\bm{\lambda}}_{i}^{H}\\!+\\!{\bm{N}}_{k},$
(10)

where ${\bm{N}}_{k}$ is the noise matrix consisting of i.i.d. Gaussian
elements with zero mean and variance $N_{0}$. Then, the MMSE estimate of
${\bm{g}}_{ik}^{(\text{d})}$ is given by

${\hat{\bm{g}}_{ik}^{(\text{d})}}=\frac{\sqrt{p_{{\text{p}},i}v_{ik}^{(\text{d})}}}{{\sum\limits_{j\in{\cal
X}_{k}}p_{{\text{p}},j}v_{jk}^{(\text{d})}+N_{0}}}{\bm{Y}}_{k}^{(\text{d})}{\bm{\lambda}}_{i},\forall
i,k\in{\cal K},$ (11)

Denote
${\bm{g}}_{ik}^{(\text{d})}={\hat{\bm{g}}_{ik}^{(\text{d})}}+{\tilde{\bm{g}}_{ik}^{(\text{d})}}$.
Then, as mentioned above, ${\hat{\bm{g}}_{ik}^{(\text{d})}}$ is statistically
independent of ${\tilde{\bm{g}}_{ik}^{(\text{d})}}$, and the distributions of
them are given by

${\hat{\bm{g}}_{ik}^{(\text{d})}}\sim{\cal
CN}({\bm{0}},\mu_{ik}^{(\text{d})}{\bm{I}}_{M}),{\tilde{\bm{g}}_{ik}^{(\text{d})}}\sim{\cal
CN}({\bm{0}},\epsilon_{ik}^{(\text{d})}{\bm{I}}_{M}),$ (12)

where

$\mu_{ik}^{(\text{d})}=\frac{p_{{\text{p}},i}v_{ik}^{(\text{d})}}{\sum\limits_{j\in{\cal
X}_{k}}p_{{\text{p}},j}v_{jk}^{(\text{d})}+N_{0}},\\\
\epsilon_{ik}^{(\text{d})}=1-\mu_{ik}^{(\text{d})}.$

Similarly, we have the MMSE estimate of ${\bm{g}}_{nk}^{(\text{c})}$ as
follows

${\hat{\bm{g}}_{nk}^{(\text{c})}}=\frac{\sqrt{q_{{\text{p}},n}v_{nk}^{(\text{c})}}}{q_{{\text{p}},n}v_{nk}^{(\text{c})}+N_{0}}{\bm{Y}}_{k}^{(\text{d})}{\bm{\omega}}_{n},\forall
n\in{\cal N},k\in{\cal K}.$ (13)

Denote the channel estimation error vector by
${\tilde{\bm{g}}_{nk}^{(\text{c})}}$, then, ${\hat{\bm{g}}_{nk}^{(\text{c})}}$
and ${\tilde{\bm{g}}_{nk}^{(\text{c})}}$ are statistically independent and
satisfy

${\hat{\bm{g}}_{nk}^{(\text{c})}}\sim{\cal
CN}({\bm{0}},\mu_{nk}^{(\text{c})}{\bm{I}}_{M}),{\tilde{\bm{g}}_{nk}^{(\text{c})}}\sim{\cal
CN}({\bm{0}},\epsilon_{nk}^{(\text{c})}{\bm{I}}_{M}),$ (14)

where

$\mu_{nk}^{(\text{c})}=\frac{q_{{\text{p}},n}v_{nk}^{(\text{c})}}{q_{{\text{p}},n}v_{nk}^{(\text{c})}+N_{0}},\epsilon_{nk}^{(\text{c})}=1-\mu_{nk}^{(\text{c})}.$
(15)

From (4) and (13), it can be observed that the estimation of
${\bm{h}}_{n}^{(\text{c})}$ and ${\bm{g}}_{nk}^{(\text{c})}$ is only affected
by pilot noise. The pilot interference from other CUs and DUs disappear
completely due to the orthogonality of the pilots. As for the estimation of
${\bm{h}}_{k}^{(\text{d})}$ and ${\bm{g}}_{ik}^{(\text{d})}$, it is clear from
(7) and (11) that apart from the effect of pilot noise, it is also affected by
pilot contamination due to PR among D2D pairs.

##  IV Achievable Rate Analysis

In this section, we analyze the achievable rate of both cellular and D2D links
under PR. Let $\bm{\beta}_{n}^{(\text{c})}$ denote the unit norm receive
filter used by BS for detecting the signal of CU $n$, and
$\bm{\beta}_{k}^{(\text{d})}$ denote the unit norm receive filter adopted by
D2D-Rx $k$ for detecting the signal from D2D-Tx $k$. Since the receive filter
can be used either to boost the desired signal power or to eliminate
interference signal, the SINR of each link critically depends on the receive
filter that is used. In this paper, we adopt PZF receivers, which use part of
degrees of freedom for signal enhancement and the remaining degrees of freedom
for interference suppression, for signal detection at both BS and D2D-Rxs.

We assume that BS uses $b_{\text{c}}$ and $b_{\text{d}}$ degrees of freedom to
cancel the interference from the nearest $b_{\text{c}}$ cellular interferers
and the nearest $b_{\text{d}}$ D2D interferers using different orthogonal
pilots. From (7), we can obtain the following relationship

${\hat{\bm{h}}_{k}^{(\text{d})}}=\sqrt{\frac{p_{{\text{p}},k}u_{k}^{(\text{d})}}{p_{{\text{p}},i}u_{i}^{(\text{d})}}}{\hat{\bm{h}}_{i}^{(\text{d})}},\forall
k\in{\cal K},i\in{\cal X}_{k}\setminus k,$ (16)

which indicates that the estimation of channels from two different D2D-Txs
using the same pilot to BS are in the same direction. As a result, the
interference from D2D interferers applying the same pilot can be eliminated
simultaneously by using one degree of freedom. Since $\tau-N$ orthogonal
pilots are reused by $K$ D2D pairs, we divide all D2D pairs into $\tau-N$ sets
with D2D pairs in each set using the same pilot for channel estimation. Then,
we know that BS is able to cancel the interference from $b_{\text{d}}$ sets of
D2D interferers. The feasible set of $(b_{\text{c}},b_{\text{d}})$ is given by

$\left\\{(b_{\text{c}},b_{\text{d}})\\!\in\\!{\mathbb{N}}\\!\times\\!{\mathbb{N}}\mid
b_{\text{c}}\\!\leq\\!N\\!-\\!1,b_{\text{d}}\\!\leq\\!\tau\\!-\\!N,b_{\text{c}}\\!+\\!b_{\text{d}}\\!\leq\\!B\\!-\\!1\right\\}.$
(17)

The PZF filter $\bm{\beta}_{n}^{(\text{c})}$ can be obtained by normalizing
the projection of channel estimation ${\hat{\bm{h}}_{n}^{(\text{c})}}$ on the
nullspace of channel estimation vectors of cancelled interferers (refer to
(54) in Appendix A). For the sake of convenience, let ${\cal
C}_{n}^{(\text{c})}$ denote the set of uncancelled CUs when detecting
$x_{n}^{(\text{c})}$, and ${\cal D}^{(\text{c})}$ denote the set of
uncancelled DUs when detecting cellular signals.

Similarly, each D2D-Rx uses $m_{\text{c}}$ and $m_{\text{d}}$ degrees of
freedom to cancel the interference from the nearest $m_{\text{c}}$ cellular
interferers and the nearest $m_{\text{d}}$ D2D interferers using different
orthogonal pilots. $(m_{\text{c}},m_{\text{d}})$ should be in the following
set

$\left\\{\\!(\\!m_{\text{c}},m_{\text{d}})\\!\in\\!{\mathbb{N}}\\!\times\\!{\mathbb{N}}\\!\mid\\!m_{\text{c}}\\!\leq\\!N,m_{\text{d}}\\!\leq\\!\tau\\!-\\!N\\!-\\!1,m_{\text{c}}\\!+\\!m_{\text{d}}\\!\leq\\!M\\!\\!-\\!\\!1\\!\right\\}\\!.$
(18)

The PZF filter $\bm{\beta}_{k}^{(\text{d})}$ can be obtained by first
projecting channel estimation ${\hat{\bm{g}}_{kk}^{(\text{d})}}$ onto the
nullspace of channel estimation vectors of cancelled interferers, and then
normalizing the projection. Let ${\cal C}_{k}^{(\text{d})}$ and ${\cal
D}_{k}^{(\text{d})}$ respectively denote the sets of uncancelled CU and DUs
when detecting $x_{k}^{(\text{d})}$.

Compared with MMSE receivers which optimally balance signal enhancement and
interference suppression, PZF receivers are suboptimal. However, we adopt PZF
receivers in this paper due to the following advantages. First, PZF receivers
take both signal boosting and interference cancellation into account, which is
similar as MMSE receivers in concept. It has been shown that the performance
of PZF receivers in terms of system throughput can approach that of MMSE
receivers and is much better than that of maximum ratio combining (MRC) or
fully zero forcing (ZF) receivers [32]. Second, the simple structure of PZF
receivers makes the analysis of the system more tractable, and allows us to
evaluate the performance of the D2D underlaid massive MIMO system in a more
explicit way. In addition, PZF receivers can be simplified as MRC receivers by
letting $(b_{\text{c}},b_{\text{d}})=(0,0)$ or
$(m_{\text{c}},m_{\text{d}})=(0,0)$, and can also be reduced to fully ZF
receivers by letting $(b_{\text{c}},b_{\text{d}})=(N-1,\tau-N)$ or
$(m_{\text{c}},m_{\text{d}})=(N,\tau-N-1)$, which can make the analysis of
this paper more general. In order to cancel the nearest interferers and obtain
PZF receivers, BS and D2D-Rxs have to know the positions of all transmitters.
To relax this requirement and make it more practical, we can obtain PZF
receivers by cancelling the interferers with the largest large-scale fading
coefficients.

###  IV-A A lower bound on achievable rate of cellular links

Since BS only has the information of estimated channel vectors (4) and (7),
which are treated as the true CSI, using PZF receiver
$\bm{\beta}_{n}^{(\text{c})}$ for detecting $x_{n}^{(\text{c})}$, we can write
the post-processing received signal associated with CU $n$ at BS as

$\displaystyle
r_{n}^{(\text{c})}=\left({\bm{\beta}}_{n}^{(\text{c})}\right)^{H}{\bm{y}}^{(\text{c})}$
$\displaystyle=\sqrt{q_{{\text{s}},n}u_{n}^{(\text{c})}}\left({\bm{\beta}}_{n}^{(\text{c})}\right)^{H}\\!{\hat{\bm{h}}_{n}^{(\text{c})}}x_{n}^{(\text{c})}$
$\displaystyle+\left(\\!{\bm{\beta}}_{n}^{(\text{c})}\\!\right)^{H}\\!\\!\left(\sum\limits_{a\in{\cal
C}_{n}^{(\text{c})}\setminus
n}\\!\\!\\!\sqrt{q_{{\text{s}},a}u_{a}^{(\text{c})}}{\hat{\bm{h}}_{a}^{(\text{c})}}x_{a}^{(\text{c})}\right.\\!+\\!\\!\sum\limits_{i\in{\cal
D}^{(\text{c})}}\\!\\!\sqrt{p_{{\text{s}},i}u_{i}^{(\text{d})}}{\hat{\bm{h}}_{i}^{(\text{d})}}x_{i}^{(\text{d})}$
$\displaystyle+\sum\limits_{a=1}^{N}\sqrt{q_{{\text{s}},a}u_{a}^{(\text{c})}}{\tilde{\bm{h}}_{a}^{(\text{c})}}x_{a}^{(\text{c})}+\left.\sum\limits_{i=1}^{K}\sqrt{p_{{\text{s}},i}u_{i}^{(\text{d})}}{\tilde{\bm{h}}_{i}^{(\text{d})}}x_{i}^{(\text{d})}+{\bm{z}}\\!\right)\\!,$
(19)

where only the first term of the second equality is the desired
information, while the other terms represent the cochannel interference,
channel estimation error and noise, respectively. Thus, the effective SINR of
cellular link $n$ can be expressed as

$\eta_{n}^{(\text{c})}=\frac{S_{n}^{(\text{c})}}{I_{n}^{(\text{c}\rightarrow\text{c})}+I_{n}^{(\text{d}\rightarrow\text{c})}+\alpha^{(\text{c})}\left\|{\bm{\beta}}_{n}^{(\text{c})}\right\|_{2}^{2}},$
(20)

where
$S_{n}^{(\text{c})}=q_{{\text{s}},n}u_{n}^{(\text{c})}\left|\left({\bm{\beta}}_{n}^{(\text{c})}\right)^{H}{\hat{\bm{h}}}_{n}^{(\text{c})}\right|^{2}$
represents the desired signal from CU $n$,
$I_{n}^{(\text{c}\rightarrow\text{c})}$ and
$I_{n}^{(\text{d}\rightarrow\text{c})}$ respectively denote the cochannel
interference from all uncancelled cellular and D2D interferers, and they are
given by

$\displaystyle I_{n}^{(\text{c}\rightarrow\text{c})}=\sum\limits_{a\in{\cal
C}_{n}^{(\text{c})}\setminus
n}q_{{\text{s}},a}u_{a}^{(\text{c})}\left|\left({\bm{\beta}}_{n}^{(\text{c})}\right)^{H}{\hat{\bm{h}}_{a}^{(\text{c})}}\right|^{2},$
$\displaystyle I_{n}^{(\text{d}\rightarrow\text{c})}=\sum\limits_{i\in{\cal
D}^{(\text{c})}}p_{{\text{s}},i}u_{i}^{(\text{d})}\left|\left({\bm{\beta}}_{n}^{(\text{c})}\right)^{H}{\hat{\bm{h}}_{i}^{(\text{d})}}\right|^{2}.$
(21)

$\alpha^{(\text{c})}$ characterizes the effect of both channel estimation
error and noise experienced by CU $n$, and can be formulated as

$\alpha^{(\text{c})}=\sum\limits_{a=1}^{N}q_{{\text{s}},a}u_{a}^{(\text{c})}\varepsilon_{a}^{(\text{c})}+\sum\limits_{i=1}^{K}p_{{\text{s}},i}u_{i}^{(\text{d})}\varepsilon_{i}^{(\text{d})}+N_{0}.$
(22)

In [9] and [22], the asymptotic uplink rate of cellular links in a massive
MIMO (or D2D underlaid massive MIMO) system has been studied. As for the
system considered in this paper, we can also obtain a similar result about the
asymptotic uplink rate of CUs as shown in the following corollary. Since
corollary 1 can be analogously verified as that in [22], we omit the proof
process for brevity.

###### Corollary 1

With fixed transmit powers at all transmitters, using linear filters for
signal detection at BS, the asymptotic uplink rate of each cellular link grows
unboundedly as $B$ goes to infinity.

In fact, the number of antennas at BS is usually finite due to multiple
practical constraints. Hence, in the following of this paper, we consider a
more practical scenario where BS is equipped with large but finite numbers of
antennas. Consider the block fading model, where all channels remain unchanged
over the coherence interval with length $T$. Then, based on (5), (8) and (20),
we can derive a lower bound on the ergodic achievable rate of cellular links
as shown in the following theorem.

###### Theorem 1

Given the SINR formula in (20), the ergodic achievable rate of cellular link
$n$ is lower bounded by

$R_{n}^{(\text{c},{\text{lb}})}=\left(1-\frac{\tau}{T}\right)\log_{2}\left(1+\eta_{n}^{(\text{c},{\text{lb}})}\right),\forall
n\in{\cal N},$ (23)

where

$\eta_{n}^{(\text{c},{\text{lb}})}=\frac{q_{{\text{s}},n}\phi_{n}^{(\text{c})}}{\sum\limits_{a=1}^{N}q_{{\text{s}},a}\varphi_{an}^{(\text{c})}+\sigma^{(\text{c})}},$
(24)

and

$\displaystyle\phi_{n}^{(\text{c})}$ $\displaystyle=$
$\displaystyle(B-b_{\text{c}}-b_{\text{d}}-1)u_{n}^{(\text{c})}\delta_{n}^{(\text{c})},$
$\displaystyle\varphi_{an}^{(\text{c})}$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}u_{a}^{(\text{c})},&a\in{\cal
C}_{n}^{(\text{c})}\setminus n\\\
u_{a}^{(\text{c})}\varepsilon_{a}^{(\text{c})},&a=n\,{\text{or}}\,a\in{\cal
N}\setminus{\cal C}_{n}^{(\text{c})}\end{array}\right.,$ (27)
$\displaystyle\sigma^{(\text{c})}$ $\displaystyle=$
$\displaystyle\sum\limits_{i=1}^{K}p_{{\text{s}},i}\varphi_{i}^{(\text{d})}+N_{0},$
$\displaystyle\varphi_{i}^{(\text{d})}$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}u_{i}^{(\text{d})},&i\in{\cal
D}^{(\text{c})}\\\ u_{i}^{(\text{d})}\varepsilon_{i}^{(\text{d})},&i\in{\cal
K}\setminus{\cal D}^{(\text{c})}\end{array}\right..$ (30)

Proof: See Appendix A. $\Box$

###### Remark 1

When considering the effect of pilot length $\tau$ on
$R_{n}^{(\text{c},{\text{lb}})}$, from (24), we can find that when
$b_{\text{d}}=0$, $\eta_{n}^{(\text{c},{\text{lb}})}$ is affected by channel
estimation error of cellular link $n$, i.e., $\varepsilon_{n}^{(\text{c})}$,
and channel estimation errors of $b_{\text{c}}$ cancelled cellular links,
i.e., $\varepsilon_{a}^{(\text{c})},a\in{\cal N}\setminus{\cal
C}_{n}^{(\text{c})}$. Since CUs use orthogonal pilots for channel estimation,
the value of $\tau$ has no effect on $\varepsilon_{n}^{(\text{c})}$ and
$\varepsilon_{a}^{(\text{c})},a\in{\cal N}\setminus{\cal C}_{n}^{(\text{c})}$.
Hence, for fixed pilot transmit power, $R_{n}^{(\text{c},{\text{lb}})}$
decreases monotonically w.r.t. $\tau$. In contrast, when $b_{\text{d}}>0$,
except the effect of $\varepsilon_{n}^{(\text{c})}$ and
$\varepsilon_{a}^{(\text{c})},a\in{\cal N}\setminus{\cal C}_{n}^{(\text{c})}$,
$\eta_{n}^{(\text{c},{\text{lb}})}$ is also influenced by the estimation
errors of channels from $b_{\text{d}}$ cancelled D2D-Txs to BS, i.e.,
$\varepsilon_{i}^{(\text{d})},i\in{\cal K}\setminus{\cal D}^{(\text{c})}$. Due
to PR, increasing $\tau$ results in smaller
$\varepsilon_{i}^{(\text{d})},i\in{\cal K}\setminus{\cal D}^{(\text{c})}$, and
thereby helps increase $\eta_{n}^{(\text{c},{\text{lb}})}$. Therefore, it
would be hard to determine the monotonicity of
$R_{n}^{(\text{c},{\text{lb}})}$ w.r.t. $\tau$.

###  IV-B A lower bound on achievable rate of D2D links

To detect $x_{k}^{(\text{d})}$, the received signal at D2D-Rx $k$ after using
PZF receiver ${\bm{\beta}}_{k}^{(\text{d})}$ can be written as

$\displaystyle
r_{k}^{(\text{d})}=\left({\bm{\beta}}_{k}^{(\text{d})}\right)^{H}{\bm{y}}_{k}^{(\text{d})}$
$\displaystyle=\sqrt{p_{{\text{s}},k}v_{kk}^{(\text{d})}}\left({\bm{\beta}}_{k}^{(\text{d})}\right)^{H}\\!{\hat{\bm{g}}_{kk}^{(\text{d})}}x_{k}^{(\text{d})}$
$\displaystyle+\left(\\!{\bm{\beta}}_{k}^{(\text{d})}\\!\right)^{H}\\!\\!\left(\\!\sum\limits_{i\in{\cal
D}_{k}^{(\text{d})}\setminus
k}\\!\\!\sqrt{p_{{\text{s}},i}v_{ik}^{(\text{d})}}{\hat{\bm{g}}_{ik}^{(\text{d})}}x_{i}^{(\text{d})}\right.\\!+\\!\sum\limits_{n\in{\cal
C}_{k}^{(\text{d})}}\\!\sqrt{q_{{\text{s}},n}v_{nk}^{(\text{c})}}{\hat{\bm{g}}_{nk}^{(\text{c})}}x_{n}^{(\text{c})}$
$\displaystyle+\sum\limits_{i=1}^{K}\sqrt{p_{{\text{s}},i}v_{ik}^{(\text{d})}}{\tilde{\bm{g}}_{ik}^{(\text{d})}}x_{i}^{(\text{d})}+\left.\sum\limits_{n=1}^{N}\sqrt{q_{{\text{s}},n}v_{nk}^{(\text{c})}}{\tilde{\bm{g}}_{nk}^{(\text{c})}}x_{n}^{(\text{c})}+{\bm{n}}_{k}\\!\right)\\!\\!.$
(31)

where only the first term of the second equality is the desired signal,
while the other terms respectively denote the cochannel interference, channel
estimation error and noise. Then, the effective SINR of D2D link $k$ is

$\eta_{k}^{(\text{d})}\\!=\\!\frac{S_{k}^{(\text{d})}}{I_{k}^{(\text{d}\rightarrow\text{d})}+I_{k}^{(\text{c}\rightarrow\text{d})}+\alpha_{k}^{(\text{d})}\left\|{\bm{\beta}}_{k}^{(\text{d})}\right\|_{2}^{2}},$
(32)

where
$S_{k}^{(\text{d})}=p_{{\text{s}},k}v_{kk}^{(\text{d})}\left|\left({\bm{\beta}}_{k}^{(\text{d})}\right)^{H}{\hat{\bm{g}}}_{kk}^{(\text{d})}\right|^{2}$
denotes the desired signal from D2D-Rx $k$.
$I_{k}^{(\text{d}\rightarrow\text{d})}$,
$I_{k}^{(\text{c}\rightarrow\text{d})}$ respectively denote D2D and cellular
cochannel interference, $\alpha_{k}^{(\text{d})}$ characterizes the effect of
both channel estimation error and noise experienced by D2D-Rx $k$, and they
are given by

$\displaystyle I_{k}^{(\text{d}\rightarrow\text{d})}=\sum\limits_{i\in{\cal
D}_{k}^{(\text{d})}\setminus
k}p_{{\text{s}},i}v_{ik}^{(\text{d})}\left|\left({\bm{\beta}}_{k}^{(\text{d})}\right)^{H}{\hat{\bm{g}}_{ik}^{(\text{d})}}\right|^{2},$
$\displaystyle I_{k}^{(\text{c}\rightarrow\text{d})}=\sum\limits_{n\in{\cal
C}_{k}^{(\text{d})}}q_{{\text{s}},n}v_{nk}^{(\text{c})}\left|\left({\bm{\beta}}_{k}^{(\text{d})}\right)^{H}{\hat{\bm{g}}_{nk}^{(\text{c})}}\right|^{2},$
$\displaystyle\alpha_{k}^{(\text{d})}=\sum\limits_{i=1}^{K}p_{{\text{s}},i}v_{ik}^{(\text{d})}\epsilon_{ik}^{(\text{d})}+\sum\limits_{n=1}^{N}q_{{\text{s}},n}v_{nk}^{(\text{c})}\epsilon_{nk}^{(\text{c})}+N_{0}.$
(33)

Similarly as the cellular uplink case, we can also derive a lower bound on
the ergodic achievable rate of D2D links.

###### Theorem 2

Given the SINR formula in (32), the ergodic achievable rate of D2D link $k$ is
lower bounded by

$R_{k}^{(\text{d},{\text{lb}})}=\left(1-\frac{\tau}{T}\right)\log_{2}\left(1+\eta_{k}^{(\text{d},{\text{lb}})}\right),\forall
k\in{\cal K},$ (34)

where

$\eta_{k}^{(\text{d},{\text{lb}})}=\frac{p_{{\text{s}},k}\phi_{k}^{(\text{d})}}{\sum\limits_{i=1}^{K}p_{{\text{s}},i}\psi_{ik}^{(\text{d})}+\sigma_{k}^{(\text{d})}},$
(35)

and

$\displaystyle\phi_{k}^{(\text{d})}$ $\displaystyle=$
$\displaystyle(M-m_{\text{c}}-m_{\text{d}}-1)v_{kk}^{(\text{d})}\mu_{kk}^{(\text{d})},$
$\displaystyle\psi_{ik}^{(\text{d})}$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}v_{ik}^{(\text{d})},\quad i\in{\cal
D}_{k}^{(\text{d})}\setminus{\cal X}_{k}&\\\
v_{ik}^{(\text{d})}\epsilon_{ik}^{(\text{d})},\quad i=k\,{\text{or}}\,{\cal
K}\setminus{\cal D}_{k}^{(\text{d})}&\\\
(M\\!-\\!m_{\text{c}}\\!-\\!m_{\text{d}}\\!-\\!1)v_{ik}^{(\text{d})}\mu_{ik}^{(\text{d})}\\!+\\!v_{ik}^{(\text{d})}\epsilon_{ik}^{(\text{d})},&i\in{\cal
X}_{k}\setminus k\end{array}\right.,$ (39)
$\displaystyle\sigma_{k}^{(\text{d})}$ $\displaystyle=$
$\displaystyle\sum\limits_{n\in{\cal
C}_{k}^{(\text{d})}}q_{{\text{s}},n}v_{nk}^{(\text{c})}+\sum\limits_{n\in{\cal
N}\setminus{\cal
C}_{k}^{(\text{d})}}q_{{\text{s}},n}v_{nk}^{(\text{c})}\epsilon_{nk}^{(\text{c})}+N_{0}.$
(40)

Proof: See Appendix B. $\Box$

###### Remark 2

From (35), it can be found that for any $(m_{\text{c}},m_{\text{d}})$ in
feasible set (18), $\eta_{k}^{(\text{d},{\text{lb}})}$ is an implicit function
of $\tau$ due to PR. Increasing $\tau$ results in more accurate channel
estimations of D2D links, and thereby helps increase
$\eta_{k}^{(\text{d},{\text{lb}})}$. However, as $\tau$ increases, the number
of symbols available for data transmission becomes smaller. In Section VI, we
show by simulation results that $R_{k}^{(\text{d},{\text{lb}})}$ first
increases and then decreases w.r.t. $\tau$.

In the following, we focus on pilot scheduling and power control design based
on (23) and (34). Since large-scale fading coefficients vary slowly, the
proposed pilot scheduling and power control algorithms can be performed
periodically at a coarser frame level granularity, which will greatly decrease
the computational complexity of BS. As a result, Theorem 1 and Theorem 2 are
helpful for the following analysis.

##  V Pilot Scheduling and Power Control

Up to now, we have investigated the cannel estimation as well as achievable
rate of both cellular and D2D links in a D2D underlaid massive MIMO system
with PR. Based on the above analysis, we focus on two problems in this
section. The first problem aims to minimize the sum MSE of channel estimation
of D2D links, and the second problem aims to maximize the sum rate of all D2D
links while guaranteeing the QoS requirements of CUs.

###  V-A Pilot Power Control and Pilot Scheduling

As mentioned in Section III, the channel estimation of cellular links is only
affected by additive noise. Therefore, we assume that each CU transmits pilot
signal with the maximum power. As for D2D links, apart from the effect of
additive noise, channel estimation is also influenced by pilot contamination.
Due to the location dispersion of D2D pairs and short-distance D2D
transmission, it should be preferred that the pilot contamination can be
greatly reduced by using an effective pilot scheduling and pilot power control
algorithm. According to (12) and (III-B), the sum MSE of channel estimation of
D2D links can be written as

$\sum\limits_{k=1}^{K}{\mathbb{E}}\left\\{\left\|{{\tilde{\bm{g}}_{kk}^{(\text{d})}}}\right\|_{2}^{2}\right\\}=\sum\limits_{k=1}^{K}M\epsilon_{kk}^{(\text{d})}.$
(41)

Since orthogonal pilots
$\\{{\bm{\omega}}_{N+1},\cdots,{\bm{\omega}}_{\tau}\\}$ are reused among $K$
D2D pairs, denote the PR pattern by ${\bm{O}}\in{\mathbb{R}}^{(\tau-N)\times
K}$ with each element in ${\bm{O}}$ being binary-valued. If D2D pair $k$ is
assigned pilot
${\bm{\omega}}_{t}\in\\{{\bm{\omega}}_{N+1},\cdots,{\bm{\omega}}_{\tau}\\}$,
we have $o_{t-N,k}=1$, otherwise, we have $o_{t-N,k}=0$. Denote the pilot
transmit power vector of DUs by
${\bm{p}}_{\text{p}}\triangleq\left[{p_{{\text{p}},1},\cdots,p_{{\text{p}},K}}\right]^{T}$.
Then, aiming at minimizing (41), we arrive at the following problem

$\displaystyle\mathop{\min}\limits_{{\bm{O}},{\bm{p}}_{\text{p}}}\quad$
$\displaystyle\sum\limits_{k=1}^{K}M\epsilon_{kk}^{(\text{d})}$ (42a) s.t.
$\displaystyle 0\leq p_{{\text{p}},k}\leq\tau P_{k},\forall k\in{\cal K},$
(42b) $\displaystyle\sum\limits_{t=N+1}^{\tau}o_{t-N,k}=1,\forall k\in{\cal
K},$ (42c)

where $P_{k}$ is the maximum data transmit power of D2D-Tx $k$.
Constraints (42c) indicate that each D2D pair can be allocated only one pilot.
Note that for simplicity, $\epsilon_{kk}^{(\text{d})}$ is formulated as a
function of ${\bm{O}}$ in an implicit way in (III-B). We can also equivalently
rewrite it in an explicit way as follows

$\epsilon_{kk}^{(\text{d})}\\!=\\!\\!\sum\limits_{t=N+1}^{\tau}\\!\\!o_{t-N,k}\\!\left(\\!1-\frac{p_{{\text{p}},k}v_{kk}^{(\text{d})}}{\sum\limits_{j=1}^{K}o_{t-N,j}p_{{\text{p}},j}v_{jk}^{(\text{d})}\\!+\\!N_{0}}\\!\right),\forall
k\in{\cal K}.$ (43)

To solve problem (42), we first give the optimal condition for
${\bm{p}}_{\text{p}}$ in the following theorem.

###### Theorem 3

There always exists $\tau\leq N+K$ such that when the optimal pilot scheduling
matrix ${\bm{O}}^{\text{opt}}$ has been determined, the optimal
${\bm{p}}_{\text{p}}^{\text{opt}}$ satisfies
$p_{{\text{p}},k}^{\text{opt}}=\tau P_{k},\forall k\in{\cal K}$.

Proof: See Appendix C. $\Box$ The above theorem indicates that with a proper
$\tau$ and the optimal ${\bm{O}}^{\text{opt}}$, constraints (42b) are always
active. We can explain Theorem 3 in an intuitive way as follows. When the
number of orthogonal pilots available for DUs is appropriate (i.e., with a
relatively low PR ratio) and these pilots are allocated to D2D pairs by using
the optimal pilot scheduling scheme (a special case is $\tau=N+K$ and all D2D
pairs use different orthogonal pilots for channel estimation), ignorable pilot
contamination would be caused due to the dispersive positions of D2D pairs.
Hence, all D2D-Txs should transmit their pilot signals in the maximum power to
increase the estimation accuracy. In contrast, with a small $\tau$ (i.e., with
a relatively high PR ratio), even using the optimal pilot scheduling scheme,
pilot contamination may still influence the estimation accuracy greatly and
solving (68) may yield $p_{{\text{p}},k}^{\text{opt}}=0$. In this case, we
need to enlarge the set of pilots for DUs to decrease PR ratio.

Based on Theorem 3, in the following, we assume that D2D-Txs always transmit
pilots in the maximum power. Then, problem (42) becomes

$\displaystyle\mathop{\min}\limits_{{\bm{O}}}\quad$
$\displaystyle\sum\limits_{k=1}^{K}M\epsilon_{kk}^{(\text{d})}$ (44a) s.t.
$\displaystyle\sum\limits_{t=N+1}^{\tau}o_{t-N,k}=1,\forall k\in{\cal K},$
(44b)

which is a mixed integer programming problem. The optimal ${\bm{O}}$ can be
obtained through exhaustive search (ES). Recalling (III-B), the number of
scalar multiplication required to compute the objective function in (42) is
${\cal O}(K)$. Thus, obtaining the optimal ${\bm{O}}$ through ES involves a
complexity of ${\cal O}(K(\tau-N)^{K})$. Due to the exponential complexity, it
will be impractical to run ES when the number of D2D pairs is large.
Therefore, we propose a low complexity pilot scheduling algorithm.

To mitigate pilot contamination in a multi-cell massive MIMO system, [33]
proposed the GCPA scheme, in which a metric is defined to indicate the
interference strength among CUs and a binary matrix is used to describe the
connections of CUs. However, to obtain the binary matrix, a suboptimal
threshold needs to be found by applying iterative grid search. In this paper,
we define a continuous-valued metric $\chi_{ik}$ to evaluate the potential
interference strength between D2D pair $i$ and $k$

$\chi_{ik}=\left\\{\begin{array}[]{ll}\\!\\!0,&\\!\\!\forall k\in{\cal
K},i=k\\\
\\!\\!\ln\left(1\\!+\\!\left(\frac{v_{ik}^{(\text{d})}}{v_{kk}^{(\text{d})}}\right)^{2}\\!+\\!\left(\frac{v_{ki}^{(\text{d})}}{v_{ii}^{(\text{d})}}\right)^{2}\right),&\\!\\!\forall
k\in{\cal K},i\in{\cal K}\setminus k\\\ \end{array}\right.\\!\\!.$ (45)

Denote $\Lambda$ as the set of D2D pairs which have been allocated pilots,
then, we summarize the pilot scheduling algorithm in Algorithm 1.

Initialization:

$\Lambda=\emptyset,{\bm{O}}={\bm{0}}$. Calculate $\chi_{ik},\forall
i,k\in{\cal K}$.

Pilot Allocation:

for $j=1,\cdots,K$ do

1: $k^{{}^{\prime}}=\arg\mathop{\max}\limits_{k\in{\cal
K}\setminus\Lambda}\sum\limits_{i\in{\cal K}}\chi_{ik}$,

2:
$t^{{}^{\prime}}=\arg\mathop{\min}\limits_{t\in\\{N+1,\cdots,\tau\\}}\sum\limits_{\left\\{i|{\bm{\lambda}}_{i}={\bm{\omega}}_{t}\right\\}}\chi_{ik^{{}^{\prime}}}$,

3: ${\bm{\lambda}}_{k^{{}^{\prime}}}={\bm{\omega}}_{t^{{}^{\prime}}}$,
$o_{t^{{}^{\prime}}-N,k^{{}^{\prime}}}=1$, $\Lambda=\Lambda\cup
k^{{}^{\prime}}$.

end for

Algorithm 1 Pilot Scheduling Algorithm (PSA)

The basic idea of the PSA algorithm is that the D2D pair experiencing larger
pilot contamination possesses a higher priority for pilot allocation. The main
steps in each iteration can be explained as follows. First, D2D pair
$k^{{}^{\prime}}\in{\cal K}\setminus\Lambda$ experiencing the largest
potential interference from other DUs is selected. Then, the pilot causing the
least interference to $k^{{}^{\prime}}$ is chosen. Finally, pilot
${\bm{\omega}}_{t^{{}^{\prime}}}$ is assigned to D2D pair $k^{{}^{\prime}}$,
i.e., $o_{t^{{}^{\prime}}-N,k^{{}^{\prime}}}=1$, and $\Lambda$ is updated by
$\Lambda=\Lambda\cup k^{{}^{\prime}}$. The algorithm will be carried out for
$K$ times until all D2D pairs are allocated with pilots.

###  V-B Data Power Control

In this subsection, we aim to maximize the sum rate of all DUs while
guaranteeing the QoS requirements of CUs. Since the exact expressions of
${\mathbb{E}}\left\\{R_{n}^{(\text{c})}\right\\}$ and
${\mathbb{E}}\left\\{R_{k}^{(\text{d})}\right\\}$ are unapproachable, we use
their lower bounds (23) and (34) for replacement. Simulation results show that
the gap between the ergodic achievable rate and its lower bound is marginal,
verifying the feasibility of the approximation.

Denote the data transmit power vectors of CUs and DUs by
${\bm{q}}_{\text{s}}\triangleq\left[{q_{{\text{s}},1},\cdots,q_{{\text{s}},N}}\right]^{T}$
and
${\bm{p}}_{\text{s}}\triangleq\left[{p_{{\text{s}},1},\cdots,p_{{\text{s}},K}}\right]^{T}$,
respectively. Then, we arrive at the following problem

$\displaystyle\mathop{\max}\limits_{{\bm{q}}_{\text{s}},{\bm{p}}_{\text{s}}}\quad$
$\displaystyle\sum\limits_{k=1}^{K}R_{k}^{(\text{d},{\text{lb}})}$ (46a) s.t.
$\displaystyle\eta_{n}^{(\text{c},{\text{lb}})}\geq\gamma_{n},\forall
n\in{\cal N},$ (46b) $\displaystyle 0\leq q_{{\text{s}},n}\leq Q_{n},\forall
n\in{\cal N},$ (46c) $\displaystyle 0\leq p_{{\text{s}},k}\leq P_{k},\forall
k\in{\cal K},$ (46d)

where $\gamma_{n}$ and $Q_{n}$ respectively denote the target SINR and
the maximum data transmit power of CU $n$. From (24) and (34), we can see that
${\bm{q}}_{\text{s}}$ and ${\bm{p}}_{\text{s}}$ are coupled in the expressions
of $\eta_{n}^{(\text{c},{\text{lb}})}$ and $R_{k}^{(\text{d},{\text{lb}})}$.
Moreover, the fractional structure of SINR expressions and the log($\cdot$)
operation make (46) non-convex. Therefore, it is difficult to obtain the
optimal solution of (46). In the following, we divide the optimization into
two consecutive parts. In the first part, we optimize ${\bm{q}}_{\text{s}}$
with ${\bm{p}}_{\text{s}}$ fixed, and vice versa for the other part.

####  V-B1 Data Power Control for Cellular Links

For given ${\bm{p}}_{\text{s}}$, problem (46) can be rewritten as

$\displaystyle\mathop{\max}\limits_{{\bm{q}}_{\text{s}}}\quad$
$\displaystyle\sum\limits_{k=1}^{K}R_{k}^{(\text{d},{\text{lb}})}$ (47a) s.t.
$\displaystyle\eta_{n}^{(\text{c},{\text{lb}})}\geq\gamma_{n},\forall
n\in{\cal N},$ (47b) $\displaystyle 0\leq q_{{\text{s}},n}\leq Q_{n},\forall
n\in{\cal N}.$ (47c)

Based on (24), we can write SINR constraints (47b) in a vector form as

${\bm{q}}_{\text{s}}\succeq{\bm{F}}{\bm{q}}_{\text{s}}+{\bm{\theta}}\triangleq{\bm{\Delta}}({\bm{q}}_{\text{s}}),$
(48)

where

$\mathbf{F}=\left(\begin{array}[]{ccc}\frac{\gamma_{1}\varphi_{11}^{(\text{c})}}{\phi_{1}^{(\text{c})}}&\ldots&\frac{\gamma_{1}\varphi_{N1}^{(\text{c})}}{\phi_{1}^{(\text{c})}}\\\
\vdots&\ddots&\vdots\\\
\frac{\gamma_{N}\varphi_{1N}^{(\text{c})}}{\phi_{N}^{(\text{c})}}&\ldots&\frac{\gamma_{N}\varphi_{NN}^{(\text{c})}}{\phi_{N}^{(\text{c})}}\end{array}\right),$
$\displaystyle{\bm{\theta}}=\left[\frac{\gamma_{1}\sigma^{(\text{c})}}{\phi_{1}^{(\text{c})}},\cdots,\frac{\gamma_{N}\sigma^{(\text{c})}}{\phi_{N}^{(\text{c})}}\right]^{T},$
$\displaystyle{\bm{\Delta}}({\bm{q}}_{\text{s}})=\left[\Delta_{1}({\bm{q}}_{\text{s}}),\cdots,\Delta_{N}({\bm{q}}_{\text{s}})\right]^{T}.$
(49)

In (49), ${\bm{\theta}}$ consists of scaled cochannel interference from D2D
interferers and additive noise, and ${\bm{\Delta}}({\bm{q}}_{\text{s}})$ can
be seen as an interference function [37].

###### Remark 3

If problem (47) is feasible, the SINR constraint of each CU would hold with
equality for the optimal ${\bm{q}}_{\text{s}}^{\text{opt}}$. Equivalently, for
the vector form (48), we have
${\bm{q}}_{\text{s}}^{\text{opt}}={\bm{F}}{\bm{q}}_{\text{s}}^{\text{opt}}+{\bm{\theta}}$.
This can be readily verified by reductio. Assume that
$\eta_{n}^{(\text{c},{\text{lb}})}>\gamma_{n}$ for
${\bm{q}}_{\text{s}}^{\text{opt}}$, then, we can further increase the
objective function by decreasing $q_{{\text{s}},n}^{\text{opt}}$ slightly.

According to [37], if the spectral radius of ${\bm{F}}$ is less than 1, the
optimal power vector has the form
${\bm{q}}_{\text{s}}^{\text{opt}}=\left({\bm{I}}_{N}-{\bm{F}}\right)^{-1}{\bm{\theta}}$
. To obtain ${\bm{q}}_{\text{s}}^{\text{opt}}$, matrix inversion and spectral
radius calculation are required, resulting in high complexity. As a result, we
obtain ${\bm{q}}_{\text{s}}^{\text{opt}}$ by using the following low
complexity iterative scheme

${\bm{q}}_{\text{s}}(l+1)={\bm{\Lambda}}({\bm{q}}_{\text{s}}(l)),$ (50)

where $l$ represents the time instant, and
${\bm{\Lambda}}({\bm{q}}_{\text{s}}(l))=\left\\{\min\left\\{Q_{1},\Delta_{1}({\bm{q}}_{\text{s}}(l))\right\\},\cdots,\min\left\\{Q_{N},\Delta_{N}({\bm{q}}_{\text{s}}(l))\right\\}\right\\}$.
We denote this data power control algorithm for cellular links by DPCC. By
proving that ${\bm{\Lambda}}({\bm{q}}_{\text{s}})$ is standard [38], we can
readily verify that the DPCC algorithm converges to the optimal solution for
any initial power vector ${\bm{q}}_{\text{s}}\succeq{\bm{0}}$. The detailed
proof is given in Appendix D.

####  V-B2 Data Power Control for D2D Links

For fixed ${\bm{q}}_{\text{s}}$, we arrive at the following problem. For
notational simplicity, the constant coefficient $1-\frac{\tau}{T}$ in (34) is
omitted.

$\displaystyle\mathop{\max}\limits_{{\bm{p}}_{\text{s}}}\quad$
$\displaystyle\sum\limits_{k=1}^{K}\log_{2}\left(1+\frac{p_{{\text{s}},k}\phi_{k}^{(\text{d})}}{\sum\limits_{i=1}^{K}p_{{\text{s}},i}\psi_{ik}^{(\text{d})}+\sigma_{k}^{(\text{d})}}\right)$
(51a) s.t.
$\displaystyle\sum\limits_{k=1}^{K}p_{{\text{s}},k}\varphi_{k}^{(\text{d})}\leq\zeta,$
(51b) $\displaystyle 0\leq p_{{\text{s}},k}\leq P_{k},\forall k\in{\cal K},$
(51c)

where $\zeta=\min\left\\{\zeta_{1},\cdots,\zeta_{N}\right\\}$, and
$\zeta_{n}=\frac{1}{\gamma_{n}}q_{{\text{s}},n}\phi_{n}^{(\text{c})}-\sum\limits_{a=1}^{N}q_{{\text{s}},a}\varphi_{an}^{(\text{c})}-N_{0}$.

It can be directly seen that problem (51) is non-convex due to the fractional
expression of $\eta_{k}^{(\text{d},{\text{lb}})}$ and the $\log(\cdot)$
operation in (51a). In order to solve (51), we transform it into the
equivalent form in (69) which admits suboptimal solutions using the well-known
WMMSE approach in [39]. The details are provided in Appendix E. We summarize
the data power control algorithm for D2D links (labeled as ‘DPCD’) in
Algorithm 9.

1:Initialize $f_{k}=\sqrt{P_{k}},w_{k}=1,\forall k\in{\cal K}$, and accuracy
$\rho_{1}$.

2:repeat

3: $w_{k}^{\prime}=w_{k},\forall k\in{\cal K},$

4:
$\nu_{k}=\frac{f_{k}\sqrt{\phi_{k}^{(\text{d})}}}{f_{k}^{2}\phi_{k}^{(\text{d})}+\sum\limits_{i=1}^{K}f_{i}^{2}\psi_{ik}^{(\text{d})}+\sigma_{k}^{(\text{d})}},\forall
k\in{\cal K},$

5:
$w_{k}=\left(1-\nu_{k}f_{k}\sqrt{\phi_{k}^{(\text{d})}}\right)^{-1},\forall
k\in{\cal K},$

6: Obtain $\lambda$ by using the bisection method,

7:
$f_{k}\\!=\\!\min\left\\{\\!\\!{\sqrt{P_{k}}},\frac{w_{k}\nu_{k}\sqrt{\phi_{k}^{(\text{d})}}}{w_{k}\nu_{k}^{2}\phi_{k}^{(\text{d})}+\sum\limits_{i=1}^{K}w_{i}\nu_{i}^{2}\psi_{ki}^{(\text{d})}+\lambda\varphi_{k}^{(\text{d})}}\\!\\!\right\\},\forall
k\in{\cal K},$

8:until $\sum\limits_{k=1}^{K}\left|\ln w_{k}-\ln
w_{k}^{\prime}\right|\leq\rho_{1},$

9:$p_{{\text{s}},k}=f_{k}^{2},\forall k\in{\cal K}.$

Algorithm 2 Data Power Control Algorithm for D2D Links (DPCD)

Based on the above analysis, we can alternatively optimize
${\bm{q}}_{\text{s}}$ and ${\bm{p}}_{\text{s}}$ by iteratively applying the
proposed DPCC and DPCD algorithms, and a suboptimal solution of problem (46)
can be obtained. Let JDPC denote the joint data power control algorithm which
iteratively carries out DPCC and DPCD until convergence. For brevity, we omit
the detailed description of the JDPC Algorithm.

###  V-C Convergence Analysis for the JDPC Algorithm

Since we propose to solve the original problem (46) using the JDPC algorithm,
which operates in an iterative mechanism, it is necessary to characterize its
convergence. In each iteration, the optimal ${\bm{q}}_{\text{s}}$ is first
obtained by using the DPCC algorithm with fixed ${\bm{p}}_{\text{s}}$. Then,
for determined ${\bm{q}}_{\text{s}}$, the DPCD algorithm outputs a suboptimal
${\bm{p}}_{\text{s}}$. As a result, the sum SE of D2D links increases in each
iteration. Noting the fact that the sum SE of D2D links is always upper
bounded, we can conclude that the JDPC algorithm converges to a suboptimal
solution of problem (46).

###  V-D Implementation and Complexity Analysis

Based on the above analysis, within a coherence interval, each transmitter
first transmits its pilot with the maximum power for channel estimation, and
then transmits its signal using the power obtained by solving (46). To
successfully implement the proposed PSA and JDPC algorithms, BS requires the
information of large-scale fading coefficients
$u_{n}^{(\text{c})},u_{k}^{(\text{d})},v_{nk}^{(\text{c})},v_{ik}^{(\text{d})},\forall
n\in{\cal N},\forall i,k\in{\cal K}$.
$u_{n}^{(\text{c})},u_{k}^{(\text{d})},\forall n\in{\cal N}$ can be directly
estimated by BS, whereas $v_{nk}^{(\text{c})},v_{ik}^{(\text{d})},\forall
n\in{\cal N},\forall i,k\in{\cal K}$ need to be estimated at all D2D-Rxs and
fed back to BS. After collecting these information, BS runs the pilot
scheduling and power control algorithms, and then sends the results to all
users. Since large-scale fading coefficients vary slowly, the proposed
algorithms can be carried out at a coarser frame level granularity.

In the following, we analyze the computational complexity of the PSA and JDPC
algorithms. Since the PSA algorithm has a complexity of ${\cal O}(K^{2})$ per
iteration and it is carried out $K$ times, the total complexity is ${\cal
O}(K^{3})$. As for the JDPC algorithm, assume that it is carried out $L_{1}$
times, and in each loop $L_{2}$ and $L_{3}$ iterations are required for the
DPCC and DPCD algorithms to converge. The major complexity of the DPCC
algorithm lies in computing ${\bm{F}}{\bm{q}}_{\text{s}}$, which involves a
complexity of ${\cal O}(N^{2})$. Hence, the total complexity of the DPCC
algorithm is ${\cal O}(L_{2}K^{2})$. Since the complexity of obtaining the
Lagrange multipliers using the bisection method for accuracy $\rho$ is ${\cal
O}(\log_{2}(1/\rho))$, the overall complexity of solving problem (51) using
the DPCD algorithm is ${\cal O}(L_{3}\log_{2}(1/\rho))$. Therefore, the JDPC
algorithm has a total complexity of ${\cal
O}(L_{1}(L_{2}K^{2}+L_{3}\log_{2}(1/\rho)))$. Simulation results show that
$L_{1}$, $L_{2}$ and $L_{3}$ are small, so the proposed algorithms involve low
complexity for efficient solutions.

##  VI Simulation Results

In this section, we present simulation results to evaluate the performance of
the D2D underlaid massive MIMO system with PR. We consider a network in a
$1000\,\text{m}\times 1000\,\text{m}$ square area and all transmitters are
located uniformly in this cell. The distance between a D2D-Tx and its
associated receiver is uniformly distributed in the range of
$[0\,\text{m},D_{\max}\,\text{m}]$. For brevity, we assume equal minimum SINR
requirement for CUs, i.e., $\gamma_{n}=\gamma,\forall n\in{\cal N}$, and equal
maximum power constraint for all transmitters, i.e., $Q_{n}=P_{k}=P,\forall
n\in{\cal N},k\in{\cal K}$. According to Theorem 3, with a relatively low PR
ratio, all D2D-Txs can transmit their pilot signals in the maximum power to
increase the estimation accuracy. Therefore, in the following, we assume that
$q_{{\text{p}},n}=p_{{\text{p}},k}=\tau P,\forall n\in{\cal N},k\in{\cal K}$.
Unless otherwise specified, the other system parameters are summarized in
Table I. All simulation results are obtained by averaging over $10^{4}$
channel realizations, and each channel realization is obtained by generating a
random user distribution as well as a random set of fading coefficients.

###  VI-A Spectral Efficiency Versus the Corresponding Lower Bound

TABLE I: Simulation Parameters Maximum data transmit power $P$ | 17 dBm
---|---
Additive noise power $N_{0}$ | -100 dBm
Path loss exponent | 3.7
Standard deviation of log-normal shadowing fading | 8 dB
Accuracy $\rho$, $\rho_{1}$ | $10^{-3}$
Figure 1: Simulated cellular sum SE and its lower bound versus the number of
BS antennas with $q_{{\text{s}},n}=p_{{\text{s}},k}=P,\forall n\in{\cal
N},k\in{\cal K}$, $N=5$, $K=20$, $M=8$, $T=50$ and $D_{\max}=100$.

In this subsection, we evaluate the tightness between the simulated SE and the
corresponding lower bound of both cellular and D2D links.

First, we compare the simulated cellular sum SE with the corresponding lower
bound under different values of pilot length and different receive filters at
BS in Fig. 1. It can be seen that cellular sum SE increases with the number of
BS antennas for all considered configurations, and fully ZF receivers (i.e.,
$(b_{\text{c}},b_{\text{d}})=(N-1,\tau-N)$) greatly outperform MRC receivers
(i.e., $(b_{\text{c}},b_{\text{d}})=(0,0)$) in terms of cellular SE. As pilot
length $\tau$ grows, cellular sum SE decreases for the MRC case, while
increases for the fully ZF case. This has been explained in Remark 1.
Moreover, Fig. 1 also shows that the lower bound of cellular sum SE closely
matches the simulation for the MRC case, and almost coincides with the
simulation for the fully ZF case.

Figure 2: Simulated sum SE of D2D links and its lower bound versus pilot
length with $q_{{\text{s}},n}=p_{{\text{s}},k}=P,\forall n\in{\cal
N},k\in{\cal K}$, $N=5$, $K=20$, $B=1024$, $T=50$ and $D_{\max}=100$.

Next, in Fig. 2, we compare the simulated sum SE of D2D links with the
corresponding lower bound under different numbers of D2D-Rx antennas and
different receive filters. Note that when D2D-Rxs adopt PZF receivers for
signal detection, as mentioned in Section IV, $(m_{\text{c}},m_{\text{d}})$
should be in feasible set (18). Therefore, in Fig. 2, we assume that
$\left(m_{\text{c}},m_{\text{d}}\right)=\left(1,\min\\{\tau-6,1\\}\right)$
when $M=4$, and
$\left(m_{\text{c}},m_{\text{d}}\right)=\left(1,\min\\{\tau-6,2\\}\right)$
when $M=8$. It can be seen from this figure that PZF receivers greatly
outperform MRC receivers in terms of D2D sum SE. When $M=4$, the lower bound
on the sum SE of DUs obtained by PZF receivers is approximate to the simulated
SE obtained by MRC receivers. The gap between the simulated SE and its lower
bound is small for all considered configurations, indicating that it is
feasible to solve the data power control problem in Section V-B based on lower
bounds (24) and (34). Fig. 2 also shows that when $\tau\leq 9$, D2D sum SE
increases with pilot length, while after that D2D sum SE decreases almost
linearly with pilot length. This is because only a few orthogonal pilots are
reused among D2D pairs for a small $\tau$. In this case, the channel
estimation is significantly influenced by pilot contamination. Therefore, the
SE of D2D links can be enhanced by increasing $\tau$. However, as $\tau$
becomes large enough, the channel estimation accuracy can be hardly improved
by further enlarging $\tau$. Counterproductively, increasing pilot length
reduces the number of symbols available for data transmission, and thereby
decreases the SE of D2D links.

In the following simulation, we assume that BS applies fully ZF receivers for
signal detection. As for D2D communication, since the number of D2D-Rx
antennas may be smaller than the number of all transmitters in the cell, we
assume that D2D-Rxs apply PZF receivers for signal detection. In addition,we
set $B=1024$ and $M=8$.

###  VI-B Performance of the PSA Algorithm

Figure 3: Performances of different pilot scheduling algorithms versus pilot
length with $N=5$, $K=20$ and $D_{\max}=100$. Figure 4: Convergence behaviors
of the inner updates in the first outer iteration of the JDPC algorithm with
$N=5$, $K=20$, $T=50$, $\gamma=5$ dB and $D_{\max}=100$.

In this subsection, we evaluate the performance of the proposed PSA algorithm.
For comparison, in Fig. 3, we plot sum MSE (41) versus pilot length for
different pilot scheduling algorithms: the proposed PSA scheme, the GCPA
algorithm proposed in [33] and the random pilot scheduling (RPS). As a
benchmark, we also depict the lower bound on the sum MSE, which is obtained
when all users apply different orthogonal pilots for channel estimation. It
can be observed from Fig. 3 that the sum MSE decreases with pilot length for
all algorithms, which is consistent with intuition. Moreover, the proposed PSA
scheme approaches the lower bound quickly as $\tau$ increases, and outperforms
the other two algorithms significantly in terms of sum MSE.

###  VI-C Performance of the JDPC Algorithm

In this subsection, we evaluate the performance of the proposed JDPC
algorithm. Denote the orthogonal training scheme (i.e., $\tau=N+K$) by ‘ OT ’.

Fig. 4 and Fig. 5 illustrate the convergence behaviors of the proposed power
control algorithm under different configurations. Specifically, the left and
right panels of Fig. 4 respectively correspond to the DPCC algorithm and the
DPCD algorithm, while Fig. 5 corresponds to the JDPC algorithm. It can be seen
from Fig. 4 that the sum SE of D2D links monotonically increases during the
iterative procedure and converges rapidly for both DPCC and DPCD algorithms.
By iteratively carrying out these two algorithms to optimize the data transmit
power of CUs and DUs, the sum SE of D2D links can be effectively increased.
Fig. 5 shows that the JDPC algorithm converges after only a few iterations
(within 3 iterations for all considered configurations). This makes the
proposed data power control algorithm suitable for practical applications.

Figure 5: Convergence behaviors of the JDPC algorithm with $N=5$, $K=20$,
$T=50$, $\gamma=5$ dB and $D_{\max}=100$. Figure 6: Sum SE of the system
versus the maximum distance between a D2D pair with
$(m_{\text{c}},m_{\text{d}})=(1,2)$, $N=5$, $K=20$, $T=50$ and $\gamma=5$ dB.

In Fig. 6, the sum SE of the system versus the maximum distance between a D2D
pair is depicted. We compare the proposed scheme with the iterative channel
allocation and power control (ICAPC) algorithm proposed in [40]. Note that
[40] considered a D2D underlaid cellular system with each transceiver equipped
with one antenna. Therefore, to be fair, we set the same system configurations
when simulating the ICAPC algorithm. We also include the sum SE curve obtained
by using classical massive MIMO communication as a benchmark in Fig. 6. To
obtain this benchmark, instead of applying direct communication between D2D
pairs, data signals of D2D-Txs are first forwarded to BS and then sent to
D2D-Rxs. Let $R_{n}^{(\text{c})}$ denote the uplink SE of cellular link $n$,
and $R_{k}^{(\text{d})}$ denote the SE of the transmission from D2D-Tx $k$ to
D2D-Rx $k$ with BS working as the relay. Then, the benchmark can be obtained
by maximizing
$R=\sum\limits_{n=1}^{N}R_{n}^{(\text{c})}+\sum\limits_{k=1}^{K}R_{k}^{(\text{d})}$.
When fully ZF receivers are applied at the BS, the maximum $R$ can be obtained
by letting all transmitters using the maximum power for signal transmission.
As expected, Fig. 6 shows that the sum SE of the system decreases with
$D_{\max}$ when D2D communication is adopted. Since the ICAPC algorithm
assumed perfect CSI and aimed to maximize the sum SE of DUs based on
instantaneous CSI, the ICAPC algorithm outperforms the proposed scheme
slightly in terms of system throughput when OT is adopted. However, with PR
among DUs, the proposed scheme can obviously increase the system SE. Moreover,
we can also conclude from Fig. 6 that due to the property of short-distance
transmission, D2D communication can help improve the sum SE of a massive MIMO
system significantly especially when $D_{\max}$ is small.

Figure 7: Sum SE of D2D links versus the number of D2D pairs under different
values of coherence block length with $(m_{\text{c}},m_{\text{d}})=(1,2)$,
$\tau=10$, $T=50$, $\gamma=5$ dB and $D_{\max}=100$. Figure 8: Sum SE of D2D
links versus the coherence block length with
$(m_{\text{c}},m_{\text{d}})=(1,2)$, $N=5$, $\gamma=5$ dB and $D_{\max}=100$.

Fig. 7 depicts the sum SE of D2D links versus the number of D2D pairs under
different values of $N$. It can be seen that the sum SE of DUs first increases
with $K$ and then approaches a saturation point when $K>150$. This is because
we set $\tau=10$, as $K$ grows, PR ratio increases, resulting in more pilot
contamination for D2D channel estimations. Moreover, since a larger $N$
results in more cochannel interference and requires longer pilot overhead, as
Fig. 7 shows, the sum SE of DUs decreases with $N$.

Figure 9: Sum SE of D2D links versus the target SINR of CUs under different
$(m_{\text{c}},m_{\text{d}})$ with $N=5$, $K=20$ , $\tau=10$, $T=50$ and
$D_{\max}=100$.

The impact of the coherence block length $T$ on the sum SE of DUs is
investigated in Fig. 8. As expected, the sum SE of DUs increases with $T$ for
all considered configurations. For both the proposed scheme with OT and the
ICAPC algorithm, as $T$ grows, the sum SE of DUs with $K=40$ is first lower
and then higher than that of the $K=20$ case. This is because with OT, in the
small $T$ regime, the sum SE of DUs is mainly influenced by $\tau$, while in
the high $T$ regime, it is mainly affected by $K$. Moreover, Fig. 8 also shows
that compared with OT, the sum SE of DUs can be greatly increased by PR
especially when $T$ is small.

Fig. 9 illustrates how the target SINR of CUs and PZF parameters
$(m_{\text{c}},m_{\text{d}})$ affect the sum SE of D2D links. Several
observations can be made from this figure. First, the sum SE of D2D links
decreases with $\gamma$, and the loss in D2D sum SE resulted from the increase
of $\gamma$ reduces as $m_{\text{c}}$ grows from 0 to 2. Second, PZF receivers
(i.e., $m_{\text{c}}+m_{\text{d}}>0$) outperform the MRC receiver (i.e.,
$m_{\text{c}}=m_{\text{d}}=0$) in terms of D2D sum SE for all different
choices of $(m_{\text{c}},m_{\text{d}})$. Moreover, for the considered
configuration, the maximum sum SE of D2D links is obtained when
$(m_{\text{c}},m_{\text{d}})=(1,2)$. Further increasing the degrees of freedom
for interference suppression results in decrease of D2D sum SE since less
degrees of freedom are left for signal enhancement.

##  VII Conclusions

In this paper, we consider a D2D underlaid massive MIMO system. Due to the
fact that D2D pairs are usually located dispersively and conduct short-
distance transmission in low power, letting several D2D pairs far from each
other use the same pilot for channel estimation would be feasible and
beneficial. Hence, we allow PR among DUs to reduce pilot overhead. Based on
this setup, we first investigate the channel estimation under PR and derive a
lower bound on the ergodic achievable rate of each link. To mitigate pilot
contamination caused by PR, we develop a pilot scheduling algorithm under the
criterion of minimizing the sum MSE of channel estimation of D2D links. In
addition, we also maximize the sum rate of all D2D links based on the large-
scale fading coefficients instead of the instantaneous CSI. An iterative power
control algorithm is proposed to obtain a suboptimal solution. Simulation
results show that the effect of pilot contamination can be decreased greatly
by exploiting the proposed pilot scheduling algorithm, and the PR scheme can
provide significant performance gains over the conventional orthogonal
training scheme in terms of system SE.

## Appendix A
Proof of Theorem 1

By definition,
$\left(\\!{\bm{\beta}}_{n}^{(\text{c})}\\!\right)^{H}{\hat{\bm{h}}}_{n}^{(\text{c})}$,
$\left(\\!{\bm{\beta}}_{n}^{(\text{c})}\\!\right)^{H}{\hat{\bm{h}}}_{a}^{(\text{c})}$
and
$\left(\\!{\bm{\beta}}_{n}^{(\text{c})}\\!\right)^{H}{\hat{\bm{h}}}_{i}^{(\text{d})}$
are respectively the projections of channel estimation vectors
${\hat{\bm{h}}}_{n}^{(\text{c})}$, ${\hat{\bm{h}}}_{a}^{(\text{c})}$ and
${\hat{\bm{h}}}_{i}^{(\text{d})}$ onto PZF receiver
${\bm{\beta}}_{n}^{(\text{c})}$. As mentioned in Section III, channel
estimation of cellular links is only affected by additive noise. Hence,
${\hat{\bm{h}}_{n}^{(\text{c})}}$ is independent of
${\hat{\bm{h}}_{a}^{(\text{c})}}$ and ${\hat{\bm{h}}_{i}^{(\text{d})}}$,
$\forall n\\!\in\\!{\cal N},a\\!\in\\!{\cal N}\setminus n,i\\!\in\\!{\cal K}$.
Then, from (20), we can find that the desired signal of CU $n$
$S_{n}^{(\text{c})}$ doesn’t dependent on cochannel cellular and D2D
interference $I_{n}^{(\text{c}\rightarrow\text{c})}$ and
$I_{n}^{(\text{d}\rightarrow\text{c})}$. Using the convexity of
$\log_{2}\left(1+\frac{1}{x}\right)$ ($\forall x>0$) and applying the Jensen’s
inequality, we have

$\displaystyle{\mathbb{E}}\left(R_{n}^{(\text{c})}\right)\geq
R_{n}^{(\text{c},{\text{lb}})}$
$\displaystyle\triangleq\left(1-\frac{\tau}{T}\right)\log_{2}\left(1+\left({\mathbb{E}}\left\\{\frac{1}{\eta_{n}^{(\text{c})}}\right\\}\right)^{-1}\right),\forall
n\\!\in\\!{\cal N},\;\;$ (52)

where $R_{n}^{(\text{c})}$ and $R_{n}^{(\text{c},{\text{lb}})}$ denote the
instantaneous rate and a lower bound on the ergodic achievable rate of
cellular link $n$, respectively. Based on (20) and the above analysis, we have

${\mathbb{E}}\left\\{\frac{1}{\eta_{n}^{(\text{c})}}\right\\}={\mathbb{E}}\left\\{\frac{1}{S_{n}^{(\text{c})}}\right\\}\left({\mathbb{E}}\left\\{I_{n}^{(\text{c}\rightarrow\text{c})}\right\\}+{\mathbb{E}}\left\\{I_{n}^{(\text{d}\rightarrow\text{c})}\right\\}+\alpha^{(\text{c})}\right).$
(53)

In order to obtain an explicit expression of (53), we need to calculate
${\mathbb{E}}\left\\{\frac{1}{S_{n}^{(\text{c})}}\right\\}$,
${\mathbb{E}}\left\\{I_{n}^{(\text{c}\rightarrow\text{c})}\right\\}$ and
${\mathbb{E}}\left\\{I_{n}^{(\text{d}\rightarrow\text{c})}\right\\}$. For
convenience, denote $\bm{G}=[\bm{G}_{1},\bm{G}_{2}]\in{\mathbb{C}}^{B\times
B}$, where the columns of $\bm{G}$ form an orthogonal basis of the $B$
dimensional space. Specifically, the columns of
$\bm{G}_{2}\in{\mathbb{C}}^{B\times(b_{\text{c}}+b_{\text{d}})}$ form an
orthogonal basis of the subspace spanned by channel estimation vectors of
cancelled interferers, and each column of
$\bm{G}_{1}\in{\mathbb{C}}^{B\times(B-b_{\text{c}}-b_{\text{d}})}$ is
orthogonal to $\bm{G}_{2}$. Then, the unit norm PZF receiver
${\bm{\beta}}_{n}^{(\text{c})}$ can be chosen as

${\bm{\beta}}_{n}^{(\text{c})}=\frac{\bm{G}_{1}\bm{G}_{1}^{H}{\hat{\bm{h}}}_{n}^{(\text{c})}}{\left\|\bm{G}_{1}\bm{G}_{1}^{H}{\hat{\bm{h}}}_{n}^{(\text{c})}\right\|_{2}}.$
(54)

By expressing ${\hat{\bm{h}}}_{n}^{(\text{c})}$ as the sum of projections of
${\hat{\bm{h}}}_{n}^{(\text{c})}$ onto $\bm{G}_{1}$ and $\bm{G}_{2}$, we have

$\displaystyle\left({\bm{\beta}}_{n}^{(\text{c})}\right)^{H}{\hat{\bm{h}}}_{n}^{(\text{c})}$
$\displaystyle=$
$\displaystyle\left({\bm{\beta}}_{n}^{(\text{c})}\right)^{H}\left(\bm{G}_{1}\bm{G}_{1}^{H}{\hat{\bm{h}}}_{n}^{(\text{c})}+\bm{G}_{2}\bm{G}_{2}^{H}{\hat{\bm{h}}}_{n}^{(\text{c})}\right)$
(55) $\displaystyle=$
$\displaystyle\left\|\bm{G}_{1}\bm{G}_{1}^{H}{\hat{\bm{h}}}_{n}^{(\text{c})}\right\|_{2}$
$\displaystyle=$
$\displaystyle\left\|\bm{G}_{1}^{H}{\hat{\bm{h}}}_{n}^{(\text{c})}\right\|_{2},$

which indicates that
$\left({\bm{\beta}}_{n}^{(\text{c})}\right)^{H}{\hat{\bm{h}}}_{n}^{(\text{c})}$
equals to the norm of the projection of ${\hat{\bm{h}}}_{n}^{(\text{c})}$ on
$\bm{G}_{1}$. Obviously, ${\hat{\bm{h}}}_{n}^{(\text{c})}$ is independent of
each column of $\bm{G}_{1}$. Hence, each element of
$\bm{G}_{1}^{H}{\hat{\bm{h}}}_{n}^{(\text{c})}$ follows i.i.d.
${\cal{CN}}(0,\delta_{n}^{(\text{c})})$. According to [41, Section 2.1.6],
$\frac{1}{\delta_{n}^{(\text{c})}}\left\|\bm{G}_{1}^{H}{\hat{\bm{h}}}_{n}^{(\text{c})}\right\|_{2}^{2}\sim{\cal
W}_{1}(B-b_{\text{c}}-b_{\text{d}},1)$ is a central complex Wishart random
variable with $B-b_{\text{c}}-b_{\text{d}}$ degrees of freedom. Then, we have

${\mathbb{E}}\left\\{\frac{1}{S_{n}^{(\text{c})}}\right\\}=\frac{1}{q_{{\text{s}},n}u_{n}^{(\text{c})}(B-b_{\text{c}}-b_{\text{d}}-1)\delta_{n}^{(\text{c})}}.$
(56)

Since ${\bm{\beta}}_{n}^{(\text{c})}$ is an unit norm vector and is
independent of ${\hat{\bm{h}}}_{a}^{(\text{c})},\forall a\in{\cal
C}_{n}^{(\text{c})}\setminus n$,
$\left({\bm{\beta}}_{n}^{(\text{c})}\right)^{H}{\hat{\bm{h}}_{a}^{(\text{c})}}$
is the linear combination of complex Gaussian random variables and follows
${\cal{CN}}(0,\delta_{a}^{(\text{c})})$. Hence,
$\left|\left({\bm{\beta}}_{n}^{(\text{c})}\right)^{H}{\hat{\bm{h}}_{a}^{(\text{c})}}\right|^{2}$
follows exponential distribution, i.e.,
$\frac{1}{\delta_{a}^{(\text{c})}}\left|\left({\bm{\beta}}_{n}^{(\text{c})}\right)^{H}{\hat{\bm{h}}_{a}^{(\text{c})}}\right|^{2}\sim{\text{Exp}}(1),\forall
a\in{\cal C}_{n}^{(\text{c})}\setminus n$. Similarly,
$\frac{1}{\delta_{i}^{(\text{d})}}\left|\left({\bm{\beta}}_{n}^{(\text{c})}\right)^{H}{\hat{\bm{h}}_{i}^{(\text{d})}}\right|^{2}\sim{\text{Exp}}(1),\forall
i\in{\cal D}^{(\text{c})}$. Then, we have

$\displaystyle{\mathbb{E}}\left\\{I_{n}^{(\text{c}\rightarrow\text{c})}\right\\}$
$\displaystyle=$ $\displaystyle\sum\limits_{a\in{\cal
C}_{n}^{(\text{c})}\setminus
n}q_{{\text{s}},a}u_{a}^{(\text{c})}{\mathbb{E}}\left\\{\left|\left({\bm{\beta}}_{n}^{(\text{c})}\right)^{H}{\hat{\bm{h}}_{a}^{(\text{c})}}\right|^{2}\right\\}$
$\displaystyle=$ $\displaystyle\sum\limits_{a\in{\cal
C}_{n}^{(\text{c})}\setminus
n}q_{{\text{s}},a}u_{a}^{(\text{c})}\delta_{a}^{(\text{c})},$
$\displaystyle{\mathbb{E}}\left\\{I_{n}^{(\text{d}\rightarrow\text{c})}\right\\}$
$\displaystyle=$ $\displaystyle\sum\limits_{i\in{\cal
D}^{(\text{c})}}p_{{\text{s}},i}u_{i}^{(\text{d})}{\mathbb{E}}\left\\{\left|\left({\bm{\beta}}_{n}^{(\text{c})}\right)^{H}{\hat{\bm{h}}_{i}^{(\text{d})}}\right|^{2}\right\\}$
(57) $\displaystyle=$ $\displaystyle\sum\limits_{i\in{\cal
D}^{(\text{c})}}p_{{\text{s}},i}u_{i}^{(\text{d})}\delta_{i}^{(\text{d})}.$

Substituting (56) and (57) into (52) yields lower bound (23).

## Appendix B
Proof of Theorem 2

Unlike the channel estimation of cellular links which is only affected by
additive noise, pilot contamination exists when estimating D2D channels due to
PR. As a result, the estimation of ${\bm{g}}_{kk}^{(\text{d})}$ is not only
affected by noise, but also influenced by pilot contamination caused by DUs
$i\in{\cal X}_{k}\setminus k$. From (11), we can obtain the following
relationship

${\hat{\bm{g}}_{ik}^{(\text{d})}}=\sqrt{\frac{\mu_{ik}^{(\text{d})}}{\mu_{kk}^{(\text{d})}}}{\hat{\bm{g}}_{kk}^{(\text{d})}},\forall
i\in{\cal X}_{k}\setminus k.$ (58)

Denote

$\displaystyle I_{k}^{(\text{d}\rightarrow\text{d})}$ $\displaystyle=$
$\displaystyle{\hat{I}}_{k}^{(\text{d}\rightarrow\text{d})}+{\tilde{I}}_{k}^{(\text{d}\rightarrow\text{d})},$
$\displaystyle{\hat{I}}_{k}^{(\text{d}\rightarrow\text{d})}$ $\displaystyle=$
$\displaystyle\sum\limits_{i\in{\cal D}_{k}^{(\text{d})}\setminus{\cal
X}_{k}}p_{{\text{s}},i}v_{ik}^{(\text{d})}\left|\left({\bm{\beta}}_{k}^{(\text{d})}\right)^{H}{\hat{\bm{g}}_{ik}^{(\text{d})}}\right|^{2},$
$\displaystyle{\tilde{I}}_{k}^{(\text{d}\rightarrow\text{d})}$
$\displaystyle=$ $\displaystyle\sum\limits_{i\in{\cal X}_{k}\setminus
k}p_{{\text{s}},i}v_{ik}^{(\text{d})}\left|\left({\bm{\beta}}_{k}^{(\text{d})}\right)^{H}{\hat{\bm{g}}_{ik}^{(\text{d})}}\right|^{2}.$
(59)

Then, from (58), we have

$\frac{{\tilde{I}}_{k}^{(\text{d}\rightarrow\text{d})}}{S_{k}^{(\text{d})}}=\frac{\sum\limits_{i\in{\cal
X}_{k}\setminus
k}p_{{\text{s}},i}v_{ik}^{(\text{d})}\mu_{ik}^{(\text{d})}}{p_{{\text{s}},k}v_{kk}^{(\text{d})}\mu_{kk}^{(\text{d})}}.$
(60)

Analogous to the proof process in Appendix A, we can easily verify that
$S_{k}^{(\text{d})}$ is independent of
${\hat{I}}_{k}^{(\text{d}\rightarrow\text{d})}$ and
$I_{k}^{(\text{c}\rightarrow\text{d})}$. Hence, the relationship between the
ergodic achievable rate of D2D link $k$ and its lower bound can be expressed
as

$\displaystyle{\mathbb{E}}\left(R_{k}^{(\text{d})}\right)\\!\geq\\!R_{k}^{(\text{d},{\text{lb}})}$
$\displaystyle\triangleq\left(1-\frac{\tau}{T}\right)\log_{2}\left(1\\!+\\!\left({\mathbb{E}}\left\\{\frac{1}{\eta_{k}^{(\text{d})}}\right\\}\right)^{-1}\right),\forall
k\\!\in\\!{\cal K},$ (61)

where

$\displaystyle{\mathbb{E}}\left\\{\frac{1}{\eta_{k}^{(\text{d})}}\right\\}$
$\displaystyle=$
$\displaystyle{\mathbb{E}}\left\\{\\!\frac{1}{S_{k}^{(\text{d})}}\\!\right\\}\\!\left({\mathbb{E}}\left\\{{\hat{I}}_{k}^{(\text{d}\rightarrow\text{d})}\right\\}+{\mathbb{E}}\left\\{I_{k}^{(\text{c}\rightarrow\text{d})}\right\\}+\alpha_{k}^{(\text{d})}\right)$
(62) $\displaystyle+$
$\displaystyle\frac{{\tilde{I}}_{k}^{(\text{d}\rightarrow\text{d})}}{S_{k}^{(\text{d})}}.$

Besides, we also have

$\displaystyle{\mathbb{E}}\left\\{\frac{1}{S_{k}^{(\text{d})}}\right\\}=\frac{1}{p_{{\text{s}},k}v_{kk}^{(\text{d})}(M-m_{\text{c}}-m_{\text{d}}-1)\mu_{kk}^{(\text{d})}},$
$\displaystyle{\mathbb{E}}\left\\{I_{k}^{(\text{c}\rightarrow\text{d})}\right\\}=\sum\limits_{n\in{\cal
C}_{k}^{(\text{d})}}q_{{\text{s}},n}v_{nk}^{(\text{c})}\mu_{nk}^{(\text{c})},$
$\displaystyle{\mathbb{E}}\left\\{{\hat{I}}_{k}^{(\text{d}\rightarrow\text{d})}\right\\}=\sum\limits_{i\in{\cal
D}_{k}^{(\text{d})}\setminus{\cal
X}_{k}}p_{{\text{s}},i}v_{ik}^{(\text{d})}\mu_{ik}^{(\text{d})}.$ (63)

Substituting (60) and (63) into (61) yields lower bound (34).

## Appendix C
Proof of Theorem 3

When the optimal pilot scheduling scheme has been determined to allocate
pilots $\\{{\bm{\omega}}_{N+1},\cdots,{\bm{\omega}}_{\tau}\\}$ to DUs, problem
(42) can be equivalently written as

$\displaystyle\mathop{\max}\limits_{{\bm{p}}_{\text{p}}}\quad$
$\displaystyle\sum\limits_{k=1}^{K}\frac{U_{k}({\bm{p}}_{\text{p}})}{V_{k}({\bm{p}}_{\text{p}})}=\sum\limits_{k=1}^{K}\frac{p_{{\text{p}},k}v_{kk}^{(\text{d})}}{\sum\limits_{i\in{\cal
X}_{k}}p_{{\text{p}},i}v_{ik}^{(\text{d})}+N_{0}}$ (64a) s.t. $\displaystyle
0\leq p_{{\text{p}},k}\leq\tau P_{k},\forall k\in{\cal K}.$ (64b)

Obviously, problem (64) is a non-convex sum-of-ratios optimization, which
aims to maximize the summation of fractional functions. The parametric
algorithm is often adopted to solve this kind of problem if the numerator of
each summation term is concave and the denominator of each summation term is
convex[42], [43]. Since $U_{k}({\bm{p}}_{\text{p}})$ and
$V_{k}({\bm{p}}_{\text{p}})$ are both affine w.r.t. ${\bm{p}}_{\text{p}}$ for
any $k$, we can optimally solve problem (64) using the parametric algorithm.
According to [42], (64) can be equivalently transformed to the following
problem

$\displaystyle\mathop{\max}\limits_{{\bm{p}}_{\text{p}},{\bm{\xi}}}\quad$
$\displaystyle\sum\limits_{k=1}^{K}\xi_{k}$ (65a) s.t.
$\displaystyle\frac{U_{k}({\bm{p}}_{\text{p}})}{V_{k}({\bm{p}}_{\text{p}})}\geq\xi_{k},\forall
k\in{\cal K},$ (65b) $\displaystyle 0\leq p_{{\text{p}},k}\leq\tau
P_{k},\forall k\in{\cal K},$ (65c)

where ${\bm{\xi}}\\!=\\!(\xi_{1},\cdots,\xi_{K})^{T}$. In order to solve
the problem in (65), we resort to alternating optimization by using the
following Lemma. Applying [42, Lemma 2.1], we obtain

###### Lemma 1

If $({\bm{p}}_{\text{p}}^{\text{opt}},{\bm{\xi}}^{\text{opt}})$ is the optimal
solution of the above maximization problem in (65), then there exists
${\bm{\kappa}}^{\text{opt}}=(\kappa_{1}^{\text{opt}},\cdots,\kappa_{K}^{\text{opt}})^{T}$
such that ${\bm{p}}_{\text{p}}^{\text{opt}}$ is the optimal solution to the
following problem

$\displaystyle\mathop{\max}\limits_{{\bm{p}}_{\text{p}}}\quad$
$\displaystyle\sum\limits_{k=1}^{K}\kappa_{k}^{\text{opt}}\left(U_{k}({\bm{p}}_{\text{p}})-\xi_{k}^{\text{opt}}V_{k}({\bm{p}}_{\text{p}})\right)$
(66a) s.t. $\displaystyle 0\leq p_{{\text{p}},k}\leq\tau P_{k},\forall
k\in{\cal K}.$ (66b)

Meanwhile, the optimal values of
$({\bm{\kappa}}^{\text{opt}},{\bm{\xi}}^{\text{opt}})$ should satisfy

$\displaystyle\kappa_{k}^{\text{opt}}=\frac{1}{V_{k}({\bm{p}}_{\text{p}}^{\text{opt}})},\forall
k\in{\cal K},$ $\displaystyle
U_{k}({\bm{p}}_{\text{p}}^{\text{opt}})-\xi_{k}^{\text{opt}}V_{k}({\bm{p}}_{\text{p}}^{\text{opt}})=0,\forall
k\in{\cal K}.$ (67)

From Lemma 1, we can obtain
$({\bm{p}}_{\text{p}}^{\text{opt}},{\bm{\xi}}^{\text{opt}},{\bm{\kappa}}^{\text{opt}})$
by iteratively carrying out the following two steps until convergence: 1)
update $({\bm{\xi}},{\bm{\kappa}})$ based on (67) for given
${\bm{p}}_{\text{p}}$; 2) update ${\bm{p}}_{\text{p}}$ for given
$({\bm{\xi}},{\bm{\kappa}})$ by solving (66). In the last iteration of the
parametric algorithm, assume that
$({\bm{\xi}}^{\text{opt}},{\bm{\kappa}}^{\text{opt}})$ has been obtained.
Then, we find the final ${\bm{p}}_{\text{p}}^{\text{opt}}$ by solving (66),
which can be further reformulated as

$\displaystyle\mathop{\max}\limits_{{\bm{p}}_{\text{p}}}\quad$
$\displaystyle\\!\\!\sum\limits_{k=1}^{K}\\!\left(\\!p_{{\text{p}},k}\\!\left(\\!\kappa_{k}^{\text{opt}}v_{kk}^{(\text{d})}\\!-\\!\sum\limits_{i\in{\cal
X}_{k}}\\!\kappa_{i}^{\text{opt}}\xi_{i}^{\text{opt}}v_{ki}^{(\text{d})}\\!\right)\\!-\\!\kappa_{k}^{\text{opt}}\xi_{k}^{\text{opt}}N_{0}\\!\right)$
(68a) s.t. $\displaystyle 0\leq p_{{\text{p}},k}\leq\tau P_{k},\forall
k\in{\cal K}.$ (68b)

Since the above problem in (68) is linear w.r.t. ${\bm{p}}_{\text{p}}$,
it is directly known that $p_{{\text{p}},k}^{\text{opt}}\\!=\\!\tau P_{k}$ if
$\kappa_{k}^{\text{opt}}v_{kk}^{(\text{d})}\\!\geq\\!\sum\limits_{i\in{\cal
X}_{k}}\kappa_{i}^{\text{opt}}\xi_{i}^{\text{opt}}v_{ki}^{(\text{d})}$;
Otherwise, $p_{{\text{p}},k}^{\text{opt}}\\!=\\!0$. Notice that we have to
guarantee $p_{{\text{p}},k}^{\text{opt}}>0$ to estimate the channel vector of
D2D link $k$. When $p_{{\text{p}},k}^{\text{opt}}=0$, we can always change the
pilot assigned to D2D pair $k$ by increasing $\tau$ so that
$\kappa_{k}^{\text{opt}}v_{kk}^{(\text{d})}\geq\sum\limits_{i\in{\cal
X}_{k}}\kappa_{i}^{\text{opt}}\xi_{i}^{\text{opt}}v_{ki}^{(\text{d})}$. The
worst case is that $\tau=N+K$ and all D2D pairs use different orthogonal
pilots for channel estimation. In this case, pilot contamination vanishes and
each D2D-Tx would transmit its pilot with the maximum power to increase the
channel estimation accuracy. Therefore, Theorem 3 is proven.

## Appendix D
Convergence Proof of of the DPCC Algorithm

As discussed in [38, Theorem 2], the iterative process in (50) converges to
the optimal solution for any initial power vector
${\bm{q}}_{\text{s}}\succeq{\bm{0}}$ if function
${\bm{\Lambda}}({\bm{q}}_{\text{s}})$ is standard. Moreover, the sufficient
condition for ${\bm{\Lambda}}({\bm{q}}_{\text{s}})$ to be standard is that
${\bm{\Delta}}({\bm{q}}_{\text{s}})$ is standard. Therefore, we only need to
prove that ${\bm{\Delta}}({\bm{q}}_{\text{s}})$ is standard. Since the
elements in ${\bm{F}}$ are nonnegative, ${\bm{q}}_{\text{s}}\succeq{\bm{0}}$
and ${\bm{\theta}}\succ{\bm{0}}$, we have

$\bullet$ _Positivity:_
${\bm{\Delta}}({\bm{q}}_{\text{s}})={\bm{F}}{\bm{q}}_{\text{s}}+{\bm{\theta}}\succ{\bm{0}};$

$\bullet$ _Monotonicity:_ If
${\bm{q}}_{\text{s}}\succeq{\bm{q}}_{\text{s}}^{{}^{\prime}}$, then
${\bm{\Delta}}({\bm{q}}_{\text{s}})-{\bm{\Delta}}({\bm{q}}_{\text{s}}^{{}^{\prime}})={\bm{F}}({\bm{q}}_{\text{s}}-{\bm{q}}_{\text{s}}^{{}^{\prime}})\succeq{\bm{0}}$,
i.e.
${\bm{\Delta}}({\bm{q}}_{\text{s}})\succeq{\bm{\Delta}}({\bm{q}}_{\text{s}}^{{}^{\prime}})$;

$\bullet$ _Scalability:_ $\forall\varsigma>1$,
$\varsigma{\bm{\Delta}}({\bm{q}}_{\text{s}})=\varsigma{\bm{F}}{\bm{q}}_{\text{s}}+\varsigma{\bm{\theta}}\succ\varsigma{\bm{F}}{\bm{q}}_{\text{s}}+{\bm{\theta}}={\bm{\Delta}}(\varsigma{\bm{q}}_{\text{s}})$.

As a result, ${\bm{\Lambda}}({\bm{q}}_{\text{s}})$ is standard and the
iterative process in (50) converges to the optimal solution for any initial
power vector ${\bm{q}}_{\text{s}}\succeq{\bm{0}}$.

## Appendix E
Solving (51) via WMMSE Algorithm

For notational convenience, we denote $f_{k}=\sqrt{p_{{\text{s}},k}},\forall
k\in{\cal K}$. Then, (51) can be equivalently transformed into the following
weighted sum-MSE minimization problem

$\displaystyle\mathop{\min}\limits_{{\bm{w}},{\bm{\nu}},{\bm{f}}}\quad$
$\displaystyle\sum\limits_{k=1}^{K}\left(w_{k}e_{k}-\ln w_{k}\right)$ (69a)
s.t.
$\displaystyle\sum\limits_{k=1}^{K}f_{k}^{2}\varphi_{k}^{(\text{d})}\leq\zeta,$
(69b) $\displaystyle 0\leq f_{k}\leq\sqrt{P_{k}},\forall k\in{\cal K},$ (69c)

where $w_{k}$ is a positive weight variable, and $e_{k}$ is the mean-
square estimation error

$e_{k}=\left(\nu_{k}f_{k}\sqrt{\phi_{k}^{(\text{d})}}-1\right)^{2}+\nu_{k}^{2}\left(\sum\limits_{i=1}^{K}f_{i}^{2}\psi_{ik}^{(\text{d})}+\sigma_{k}^{(\text{d})}\right).$
(70)

To prove the equivalence between (51) and (69), we derive the optimal
${\bm{\nu}}$ and ${\bm{w}}$ by checking the first-order optimality condition
of problem (69)

$\displaystyle\nu_{k}^{\text{opt}}$ $\displaystyle=$
$\displaystyle\frac{f_{k}\sqrt{\phi_{k}^{(\text{d})}}}{f_{k}^{2}\phi_{k}^{(\text{d})}+\sum\limits_{i=1}^{K}f_{i}^{2}\psi_{ik}^{(\text{d})}+\sigma_{k}^{(\text{d})}},\forall
k\in{\cal K},$ $\displaystyle w_{k}^{\text{opt}}$ $\displaystyle=$
$\displaystyle\frac{1}{e_{k}},\forall k\in{\cal K}.$ (71)

Plugging (71) in (70) and simplifying (69), we have

$\displaystyle\mathop{\max}\limits_{\bm{f}}\quad$
$\displaystyle\sum\limits_{k=1}^{K}\log_{2}\left(1+\frac{f_{k}^{2}\phi_{k}^{(\text{d})}}{\sum\limits_{i=1}^{K}f_{i}^{2}\psi_{ik}^{(\text{d})}+\sigma_{k}^{(\text{d})}}\right)$
(72a) s.t.
$\displaystyle\sum\limits_{k=1}^{K}f_{k}^{2}\varphi_{k}^{(\text{d})}\leq\zeta,$
(72b) $\displaystyle 0\leq f_{k}\leq\sqrt{P_{k}},\forall k\in{\cal K},$ (72c)

which is equivalent to (51). The equivalence implies that a suboptimal
solution of problem (51) can be obtained by solving (69), which is easier to
handle since the objective function is convex w.r.t. each variable when the
other variables are fixed. The optimal ${\bm{w}}$ (or ${\bm{\nu}}$) can be
obtained based on (71) when ${\bm{\nu}}$ and ${\bm{f}}$ (${\bm{w}}$ and
${\bm{f}}$) are fixed. For given ${\bm{w}}$ and ${\bm{\nu}}$, to get the
optimal ${\bm{f}}$, we attach Lagrange multiplier $\lambda$ to the first
constraint of (69) and obtain the Lagrange function as follows

${\cal L}({\bm{f}},\lambda)\triangleq\sum\limits_{k=1}^{K}\left(w_{k}e_{k}-\ln
w_{k}\right)+\lambda\left(\sum\limits_{k=1}^{K}f_{k}^{2}\varphi_{k}^{(\text{d})}-\zeta\right).$
(73)

From [44], the first-order optimality condition of ${\cal L}$ w.r.t. $f_{k}$
yields

$f_{k}(\lambda)=\min\left\\{{\sqrt{P_{k}}},\frac{w_{k}\nu_{k}\sqrt{\phi_{k}^{(\text{d})}}}{w_{k}\nu_{k}^{2}\phi_{k}^{(\text{d})}+\sum\limits_{i=1}^{K}w_{i}\nu_{i}^{2}\psi_{ki}^{(\text{d})}+\lambda\varphi_{k}^{(\text{d})}}\right\\}.$
(74)

According to the complementary slackness condition, if
$\lambda_{\text{opt}}=0$ yields
$\sum\limits_{k=1}^{K}f_{k}^{2}\varphi_{k}^{(\text{d})}<\zeta$, then,
$f_{k}^{\text{opt}}=f_{k}(0)$. Otherwise, we have $\lambda_{\text{opt}}>0$ and
$\sum\limits_{k=1}^{K}f_{k}^{2}\varphi_{k}^{(\text{d})}=\zeta$. Since for any
$k\in{\cal K}$, $f_{k}(\lambda)$ strictly decreases with $\lambda$, we can
obtain $\lambda_{\text{opt}}$ using the bisection method. Then, substituting
$\lambda_{\text{opt}}$ into (74), we get $f_{k}^{\text{opt}}$.

## References

  * [1] FCC Spectrum Policy Task Force, “Report of the spectrum efficiency working group,” Nov. 2002.
  * [2] T. L. Marzetta, “Noncooperative cellular wireless with unlimited numbers of base station antennas,” _IEEE Trans. Wireless Commun._ , vol. 9, no. 11, pp. 3590–3600, Nov. 2010.
  * [3] F. Rusek, D. Persson, B. K. Lau, E. G. Larsson, T. L. Marzetta, O. Edfors, and F. Tufvesson, “Scaling up MIMO: Opportunities and challenges with very large arrays,” _IEEE Signal Process. Mag._ , vol. 30, no. 1, pp. 40–60, Jan. 2013.
  * [4] E. G. Larsson, O. Edfors, F. Tufvesson, and T. L. Marzetta, “Massive MIMO for next generation wireless systems,” _IEEE Commun. Mag._ , vol. 52, no. 2, pp. 186–195, Feb. 2014.
  * [5] H. Xie, B. Wang, F. Gao, and S. Jin, “A full-space spectrum-sharing strategy for massive MIMO cognitive radio systems,” _IEEE J. Sel. Areas Commun._ , vol. 34, no. 10, pp. 2537–2549, Oct. 2016.
  * [6] J. Liu, N. Kato, J. Ma, and N. Kadowaki, “Device-to-device communication in LTE-advanced networks: A survey,” _IEEE Commun. Surveys Tuts._ , vol. 17, no. 4, pp. 1923–1940, Dec. 2015.
  * [7] A. Asadi, Q. Wang, and V. Mancuso, “A survey on device-to-device communication in cellular networks,” _IEEE Commun. Surveys Tuts._ , vol. 16, no. 4, pp. 1801–1819, Apr. 2014.
  * [8] X. Lin, J. Andrews, A. Ghosh, and R. Ratasuk, “An overview of 3GPP device-to-device proximity services,” _IEEE Commun. Mag._ , vol. 52, no. 4, pp. 40–48, Mar. 2014.
  * [9] H. Q. Ngo, E. G. Larsson, and T. L. Marzetta, “Energy and spectral efficiency of very large multiuser MIMO systems,” _IEEE Trans. Commun._ , vol. 61, no. 4, pp. 1436–1449, Apr. 2013.
  * [10] H. Yang and T. L. Marzetta, “Performance of conjugate and zero-forcing beamforming in large-scale antenna systems,” _IEEE J. Sel. Areas Commun._ , vol. 31, no. 2, pp. 172–179, Feb. 2013.
  * [11] J. Hoydis, S. Ten Brink, and M. Debbah, “Massive MIMO in the UL/DL of cellular networks: How many antennas do we need?” _IEEE J. Sel. Areas Commun._ , vol. 31, no. 2, pp. 160–171, Feb. 2013.
  * [12] N. Liang, W. Zhang, and C. Shen, “An uplink interference analysis for massive MIMO systems with MRC and ZF receivers,” in _Proc. IEEE Wireless Commun. Netw. Conf. (WCNC)_ , 2015, pp. 310–315.
  * [13] T. Bai and R. W. Heath Jr, “Asymptotic coverage probability and rate in massive MIMO networks,” _Proc. IEEE GlobalSIP_ , pp. 1–5, 2014.
  * [14] D. Zhu, J. Wang, A. L. Swindlehurst, and C. Zhao, “Downlink resource reuse for device-to-device communications underlaying cellular networks,” _IEEE Sig. Process. Lett._ , vol. 21, no. 5, pp. 531–534, May 2014.
  * [15] S. Maghsudi and S. Stanczak, “Joint channel allocation and power control for underlay D2D transmission,” in _Proc. IEEE Int. Conf. Commun. (ICC)_ , London, UK., 2015, pp. 2091–2096.
  * [16] D. Feng, L. Lu, Y. Yuan-Wu, G. Y. Li, G. Feng, and S. Li, “Device-to-device communications underlaying cellular networks,” _IEEE Trans. Commun._ , vol. 61, no. 8, pp. 3541–3551, Aug. 2013.
  * [17] W. Zhao and S. Wang, “Resource allocation for device-to-device communication underlaying cellular networks: An alternating optimization method,” _IEEE Commun. Lett._ , vol. 19, no. 8, pp. 1398–1401, Aug. 2015.
  * [18] Y. Jiang, Q. Liu, F. Zheng, X. Gao, and X. You, “Energy efficient joint resource allocation and power control for D2D communications,” _IEEE Trans. Veh. Tech._ , vol. 65, no. 8, pp. 6119 – 6127, Aug. 2016.
  * [19] F. Wang, C. Xu, L. Song, and Z. Han, “Energy-efficient resource allocation for device-to-device underlay communication,” _IEEE Trans. Wireless Commun._ , vol. 14, no. 4, pp. 2082–2092, Apr. 2015.
  * [20] T. D. Hoang, L. B. Le, and T. Le-Ngoc, “Energy-efficient resource allocation for D2D communications in cellular networks,” in _Proc. IEEE Int. Conf. Commun._ , London, UK., 2015, pp. 2251–2256.
  * [21] X. Lin, R. W. Heath, and J. G. Andrews, “Spectral efficiency of massive MIMO systems with D2D underlay,” in _Proc. IEEE Int. Conf. Commun. (ICC)_ , London, UK., 2015, pp. 4345–4350.
  * [22] X. Lin, R. Heath, and J. Andrews, “The interplay between massive MIMO and underlaid D2D networking,” _IEEE Trans. Wireless Commun._ , vol. 14, no. 6, pp. 3337–3351, June 2015.
  * [23] J. Jose, A. Ashikhmin, T. L. Marzetta, and S. Vishwanath, “Pilot contamination and precoding in multi-cell TDD systems,” _IEEE Trans. Wireless Commun._ , vol. 10, no. 8, pp. 2640–2651, Aug. 2011.
  * [24] H. Yin, D. Gesbert, M. Filippou, and Y. Liu, “A coordinated approach to channel estimation in large-scale multiple-antenna systems,” _IEEE J. Sel. Areas Commun._ , vol. 31, no. 2, pp. 264–273, Feb. 2013.
  * [25] A. Ashikhmin and T. Marzetta, “Pilot contamination precoding in multi-cell large scale antenna systems,” in _Proc. IEEE ISIT_ , Cambridge, MA, USA, 2012, pp. 1137–1141.
  * [26] F. Fernandes, A. Ashikhmin, and T. L. Marzetta, “Inter-cell interference in noncooperative TDD large scale antenna systems,” _IEEE J. Sel. Areas Commun._ , vol. 31, no. 2, pp. 192–201, Feb. 2013.
  * [27] H. Q. Ngo and E. G. Larsson, “Evd-based channel estimation in multicell multiuser MIMO systems with very large antenna arrays,” in _Proc. IEEE ICASSP_ , Kyoto, Japan, 2012, pp. 3249–3252.
  * [28] L. You, X. Gao, X.-G. Xia, N. Ma, and Y. Peng, “Pilot reuse for massive MIMO transmission over spatially correlated rayleigh fading channels,” _IEEE Trans. Wireless Commun._ , vol. 14, no. 6, pp. 3352–3366, Feb. 2015.
  * [29] X. Liu, Q. He, Y. Li, L. Xiao, and J. Wang, “Pilot reuse for device-to-device underlay massive MIMO systems,” in _IEEE VTC-Fall_ , Sep. 2015, pp. 1–5.
  * [30] X. Liu, Y. Li, X. Li, L. Xiao, and J. Wang, “Pilot reuse and interference-aided MMSE detection for D2D underlay massive MIMO,” _IEEE Trans. Veh. Tech._ , vol. 66, no. 4, pp. 3116 – 3130, Apr. 2017.
  * [31] H. Xu, N. Huang, Z. Yang, J. Shi, B. Wu, and M. Chen, “Pilot allocation and power control in D2D underlay massive MIMO systems,” _IEEE Commun. Lett._ , vol. 21, no. 1, pp. 112–115, Jan. 2017.
  * [32] N. Jindal, J. G. Andrews, and S. Weber, “Multi-antenna communication in ad hoc networks:achieving MIMO gains with SIMO transmission,” _IEEE Trans. Commun._ , vol. 59, no. 2, pp. 529–540, Feb. 2011.
  * [33] X. Zhu, L. Dai, and Z. Wang, “Graph coloring based pilot allocation to mitigate pilot contamination for multi-cell massive MIMO systems,” _IEEE Commun. Lett._ , vol. 19, no. 10, pp. 1842–1845, Aug. 2015.
  * [34] M. Li, S. Jin, and X. Gao, “Spatial orthogonality-based pilot reuse for multi-cell massive MIMO transmission,” in _Proc. IEEE Int. Conf. WCSP_ , Oct. 2013, pp. 1–6.
  * [35] R. Mochaourab, E. Björnson, and M. Bengtsson, “Adaptive pilot clustering in heterogeneous massive mimo networks,” _IEEE Trans. Wireless Commun._ , vol. 15, no. 8, pp. 5555–5568, Aug. 2016.
  * [36] T. Kailath, A. H. Sayed, and B. Hassibi, _Linear Estimation_. Prentice Hall Upper Saddle River, NJ, 2000, vol. 1.
  * [37] R. Chen, J. G. Andrews, R. W. Heath, and A. Ghosh, “Uplink power control in multi-cell spatial multiplexing wireless systems,” _IEEE Trans. Wireless Commun._ , vol. 6, no. 7, pp. 2700–2711, July 2007.
  * [38] R. D. Yates, “A framework for uplink power control in cellular radio systems,” _IEEE J. Sel. Areas Commun._ , vol. 13, no. 7, pp. 1341–1347, Sep. 1995.
  * [39] 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, Sep. 2011.
  * [40] H. Xu, Z. Yang, N. Huang, J.-Y. Wang, J. Shi, and M. Chen, “Channel allocation and power control in D2D uplink underlaid cellular networks,” in _Proc. IEEE Globecom Workshops (GC Wkshps)_ , Washington, DC, 2016, pp. 1–6.
  * [41] A. M. Tulino and S. Verdú, _Random Matrix Theory and Wireless Communications_. Now Publishers Inc, 2004, vol. 1.
  * [42] Y. Jong, “An efficient global optimization algorithm for nonlinear sum-of-ratios problem,” _online: www. optimizationonline. org_ , 2012.
  * [43] L. Xu, G. Yu, and Y. Jiang, “Energy-efficient resource allocation in single-cell OFDM systems: Multi-objective approach,” _IEEE Trans. Wireless Commun._ , vol. 14, no. 10, pp. 5848–5858, Oct. 2015.
  * [44] S. Boyd and L. Vandenberghe, _Convex Optimization_. Cambridge University Press, 2004.