Pilot Optimization and Channel Estimation For Multiuser Massive MIMO Systems
Pilot Optimization and Channel Estimation For Multiuser Massive MIMO Systems
Pilot Optimization and Channel Estimation For Multiuser Massive MIMO Systems
i=1
h
i
x
H
it
+n
t
(1)
where h
i
C
M1
is the channel between the ith MS and
the BS, x
H
it
is the pilot symbol of the ith MS at period t and
n
t
C
M1
is the received noise vector at time slot t.
If we utilize N symbol periods for pilot transmission, the
received signal can be expressed in the following matrix form
Y = HX
H
+N (2)
where x
j
= [x
j1
, x
j2
, , x
jN
]
T
, X = [x
1
, x
2
, , x
K
],
H = [h
1
, h
2
, , h
K
] and N = [n
1
, , n
N
].
The channel is modeled as h
k
= [h
k1
, h
k2
, , h
kM
]
T
with h
km
=
g
km
h
km
, m, where
g
km
is the distant
dependent path loss and
h
km
is the fast fading part. It is
assumed that
h
km
is i.i.d random variable with
h
km
CN(0, 1), k, m, E{N
(i,j)
H
(k,m)
} = 0, i, j, k, m and each
element of N is also i.i.d zero mean circularly symmetric
complex Gaussian (ZMCSCG) random variable all with unit
variance. Furthermore, as the path loss depends on the distance
between each MS and the BS, it is assumed that
g
km
=
g
k
, m, and
g
k
, k are known a priory.
III. CONVENTIONAL PILOT ASSIGNMENT AND MMSE
CHANNEL ESTIMATION METHODS
In this section, we discuss the conventional pilot assign-
ment and MMSE channel estimation algorithm of [6].
When N = K, the overall pilot symbol of all MSs (X)
is designed in such a way that X
H
X = I
K
. The BS, then,
decouples the channels of each MS by multiplying the overall
received signal Y by X, i.e.,
Z YX = H+NX
z
k
= h
k
+ n
k
, k = 1, , K (3)
where n
k
= Nx
k
and x
k
is the kth column of X.
Now by introducing an MMSE receiver W
H
k
for the kth
MS, the estimated channel of the kth MS (
h
k
) is expressed as
h
k
= W
H
k
z
k
, k. (4)
And W
k
is designed such that the mean square error (MSE)
between
h
k
and h
k
is minimized. The MSE of the kth MS is
given by
k
=E{|
h
k
h
k
|
2
}
=E{|W
H
k
(h
k
+ n
k
) h
k
|
2
}
=E{|(W
H
k
I
M
)h
k
+ n
k
}
=g
k
(W
H
k
I
M
)(W
H
k
I
M
)
H
+
2
I
M
. (5)
The optimal W
k
is given as
k
W
H
k
= 0 W
k
=
g
k
g
k
+
2
I.
When K = 2N (i.e., N < K scenario), the overall pilot
symbols X can be designed as X = [U, U], where U
C
NN
with U
H
U = I
N
. With this settings
2
, like in the above
case, the BS will right multiply the received signal Y by X
to get
Z YX = H+NX
z
k
= h
k
+h
k+N
+ n
k
, k = 1, , N
= h
k
+h
Kk
+ n
k
, k = N + 1, , K. (6)
As we can see from this expression, the channel vector of the
kth MS is coupled with the channel vector of the (N + k)
or (K k)th MS. According to [5], this phenomena is called
pilot contamination. When there is pilot contamination, the
performance of the conventional MMSE channel estimation
degrades drastically.
2
Note that one can apply this pilot symbol design approach for any other
settings of K and N.
IV. PROPOSED PILOT OPTIMIZATION AND CHANNEL
ESTIMATION ALGORITHM
In this section, we discuss our pilot optimization and
channel estimation algorithm. To this end, we introduce the
variables u
k
and W
k
, and express the estimated channel of
the kth MS as
h
k
= W
H
k
Yu
k
, k (7)
where W
k
C
MM
and u
k
C
N1
are the variables
that will be determined in the sequel. Under these introduced
variables, the MSE between
h
k
and h
k
is given by
k
=E{|
h
k
h
k
|
2
}
=E{|W
H
k
(
K
i=1
h
i
x
H
i
+N)u
k
h
k
|
2
}
=u
H
k
(
K
i=1
g
i
x
i
x
H
i
+
2
I
N
)u
k
(W
H
k
W
k
) +g
k
I
M
(g
k
x
H
k
u
k
)W
H
k
(g
k
u
H
k
x
k
)W
k
.
It follows
k
=tr{
k
}
=u
H
k
(
K
i=1
g
i
x
i
x
H
i
+
2
I
N
)u
k
tr{(W
H
k
W
k
)} +g
k
I
M
(g
k
x
H
k
u
k
)tr{W
H
k
} (g
k
u
H
k
x
k
)tr{W
k
}.
As we can see from this expression, the MSE of each MS
depends on its path loss. And in practice, we would like the
channels of all MSs to be estimated reliably almost with the
same accuracy. Due to this reason, we consider a WSMSE
problem as our objective function, where the MSE weight of
each MS is set to the inverse of its path loss. This optimization
problem is mathematically formulated as
min
xk,uk,Wk
K
k=1
1
g
k
k
s.t x
H
k
x
k
P
k
(8)
where P
k
is the maximum transmission power at each MS.
For xed u
k
and x
k
, k, we can optimize W
k
by the
MMSE method as
W
k
=
g
k
x
H
k
u
k
K
i=1
g
i
x
H
i
u
k
u
H
k
x
i
+
2
u
H
k
u
k
I
M
. (9)
Substituting this W
k
back to
k
, we get the minimum MSE
k
as
k
=M
g
k
u
H
k
(g
2
k
x
k
x
H
k
)u
k
u
H
k
(
K
i=1
g
i
x
i
x
H
i
+
2
I
N
)u
k
. (10)
From this expression, one can notice that u
k
can be optimized
independently by solving the following problem
min
uk
k
= max
uk
u
H
k
(g
2
k
x
k
x
H
k
)u
k
u
H
k
(
K
i=1
g
i
x
i
x
H
i
+
2
I
N
)u
k
, k. (11)
Now by dening u
k
A
1/2
u
k
with A =
K
i=1
g
i
x
i
x
H
i
+
2
I, we can reformulate the above problem as
max
uk
u
H
k
A
1/2
(g
2
k
x
k
x
H
k
)A
1/2
u
k
u
H
k
u
k
.
This problem is the well known Rayleigh quotient optimiza-
tion problem and its optimal solution is given by [8], [9]
u
k
= A
1/2
g
k
x
k
u
k
= A
1/2
u
k
= A
1
g
k
x
k
. (12)
Plugging this u
k
in (10) yields
k
=Mg
k
Mg
2
k
x
H
k
A
1
x
k
. (13)
By dening V
k
x
k
x
H
k
,
H
k
g
k
I
N
, V =
blkdiag(V
1
, , V
K
) and
H = [
H
1
,
H
K
], we can ex-
press
K
k=1
1
gk
k
as
w
=
K
k=1
1
g
k
k
=
K
k=1
M(1 g
k
x
H
k
A
1
x
k
)
=
K
k=1
M Mtr{(
HV
H
H
+
2
I)
1
HV
H
H
}.
If we apply matrix inversion lemma [10], we can rewrite
(
HV
H
H
+
2
I)
1
HV
H
H
as
(
HV
H
H
+
2
I)
1
HV
H
H
= I
2
(
HV
H
H
+
2
I)
1
.
Therefore, problem (8) can be equivalently formulated as
min
Vk
tr{(
K
k=1
g
k
V
k
+
2
I
N
)
1
}
s.t tr{V
k
} P
k
,
V
k
0, Rank(V
k
) = 1, k (14)
where (.) 0 denotes a PSD constraint. When N = 1 the
rank constraint is always satised. And this solution is global
optimal solution. Furthermore, when N = K, the optimal
solution of the above problem can be obtained in closed form
by applying the well known Majorization theory of [11] and
is given by the conventional orthogonal pilot.
Now how can we solve this problem when 1 < N < K
(i.e., the most practically relevant scenario). If we relax the
rank constraint of this problem, we will get a standard convex
SDP problem which can be solved efciently with interior
point methods with polynomial time complexity [8]. However,
the drawback of this convex reformulation is that the optimal
V
k
is always a scaled diagonal matrix which is full rank.
And we are not aware of any method to get the desired rank
1 solution from the solution of the above problem. Therefore,
the critical issue will be how to get a rank one solution for
the aforementioned settings of N and K. In the following,
we provide simple iterative algorithm that will give rank one
solution for the above problem.
By employing matrix inversion lemma, we can re-express
(
K
k=1
g
k
V
k
+
2
I
N
)
1
as [10]
(
K
i=1
g
i
V
i
+
2
I
N
)
1
=(Q
k
+g
k
x
k
x
H
k
)
1
=Q
1
k
g
k
Q
1
k
x
k
x
H
k
Q
1
k
1 +g
k
x
H
k
Q
1
k
x
k
tr{(
K
i=1
g
i
V
i
+
2
I
N
)
1
} =tr{Q
1
k
}
g
k
x
H
k
Q
2
k
x
k
1 +g
k
x
H
k
Q
1
k
x
k
where Q
k
=
K
i=1,i=k
g
i
x
i
x
H
i
+
2
I
N
.
Keeping Q
k
constant, x
k
can be optimized by solving the
following problem
max
xk
g
k
x
H
k
Q
2
k
x
k
1 +g
k
x
H
k
Q
1
k
x
k
, s.t x
H
k
x
k
P
k
. (15)
Since the objective function of this problem increases as x
H
k
x
k
increases, at optimality x
H
k
x
k
= P
k
is satised. This problem
can thus be reformulated as
max
xk
g
k
x
H
k
Q
2
k
x
k
x
H
k
(
1
Pk
I
N
+g
k
Q
1
k
)x
k
. (16)
The optimal solution of this problem can be obtained like in
problem (11) and is given as x
k
=
k
F
1/2
k
x, where x is the
eigenvector corresponding to the maximum eigenvalue of the
matrix g
k
F
1/2
k
Q
2
k
F
1/2
k
with F
k
= g
k
Q
1
k
+
1
Pk
I
N
and
k
is selected such that x
H
k
x
k
= P
k
is ensured.
The proposed pilot symbol design and channel estimation
algorithm is summarized in Algorithm I.
Algorithm I
Input parameters: Set N, K,
2
and g
k
, k.
Initialize n = 0 and x
n
k
, k.
Repeat
Using x
n
k
, k, solve (16) for k = 1, , K sequentially
and update x
n
k
, k by the new x
k
, k.
Increase n by 1 (i.e., n=n+1).
Until convergence
Compute u
k
and W
k
with (12) and (9), respectively.
Using these u
k
and W
k
, compute the estimated channel
h
k
with (7).
Convergence analysis: At each iteration, as the objec-
tive function of (16) (i.e., equivalent to (14)) is non-
decreasing, this iterative algorithm is always convergent.
Computational complexity: As will be clear later in the
simulation section, this algorithm converges in n < 2
iterations. Thus, the main computational cost of this algo-
rithm arises from solving the Rayleigh quotient problem
(16) which requires the computation of the eigenvalue
decomposition of an N N matrix. And in practice as
N << K << M, the complexity of this algorithm is
almost the same as the conventional pilot symbol and
channel estimation algorithm discussed in Section III.
When K = N, it can be easily seen that this iterative
algorithm will give the optimal solution of the SDP problem
(14) when the initial x
k
, k satises X
H
X = I
K
.
We would like to point out that the above pilot optimization
and channel estimation algorithm can be applied for any M,
N and K. However, as discussed in [5] the effect of pilot
contamination (i.e., pilot reuse) is worse for massive MIMO
systems. For this reason, we believe that the current algorithm
is more relevant for massive MIMO systems than that of the
conventional MIMO systems.
V. SIMULATION RESULTS
In this section, we provide simulation results. All of our
results are obtained by averaging 20000 channel realizations
and the WSMSEs are the normalized WSMSE (i.e., achieved
WSMSE of (8) divided by KM). For this simulation, we
consider a downlink multiuser massive MIMO system with
M = 128, K = 32, P
k
= 1mw, k, the path losses are
normalized to 1 (i.e., 0 < g
k
< 1, k) and the signal to noise
ratio (SNR) is dened as
Pav
2
, where P
av
=
1
K
K
k=1
P
k
. We
also compare the performance of the proposed algorithm (i.e.,
Section IV) with that of the existing algorithm discussed in
Section III.
Unless stated otherwise, we set the overall pilot symbol
X as the rst K columns of the matrix
X for the existing
algorithm, and we use this X as an initialized pilot symbols
for the proposed algorithm, where
X = [U U] (i.e., with
pilot reuse) and Uis the normalized discrete Fourier transform
(DFT) matrix of size N.
G =
. (17)
A. Effect of SNR
In this simulation, we examine the effect of SNR on
the performance of the proposed and existing algorithms.
To this end, we set N = 16 and g
k
, k are taken from
a uniform random variable in between 0 and 1 which are
given in (17) (i.e., vectorization of G). With these settings,
we plot the achieved normalized WSMSE as a function of
SNR. As can be seen from Fig. 1, the performance of the
proposed algorithm outperforms the existing algorithm for all
SNR regions. Furthermore, for both algorithms, the normalized
WSMSE decreases as the SNR increases which is expected.
B. Effect of the number of pilot symbols N
In this simulation, we examine the effect of N on the
performance of the proposed and existing algorithms. Fig.
2 shows the normalized WSMSE of the proposed and the
existing algorithm for different settings of N and SNR when
g
k
, k are as in (17) (i.e., vectorization of G). This gure also
shows that the proposed algorithm achieves less normalized
WSMSE compared to that of the existing one for all SNR
regions when N < K. And as expected, both algorithms
achieve the same normalized WSMSE when N = K, and
their WSMSE decreases as N increases.
0 2 4 6 8 10 12 14 15
0.5
0.55
0.6
0.65
0.7
0.75
0.8
0.85
SNR (dB)
N
o
r
m
a
l
i
z
e
d
W
S
M
S
E
Existing algorithm
Proposed algorithm
Fig. 1. Comparison of the proposed algorithm (Algorithm I) and the existing
algorithm for K = 32 and N = 16
16 18 20 22 24 26 28 30 32
0.1
0.2
0.3
0.4
0.5
0.6
0.7
Number of pilot symbols (N)
N
o
r
m
a
l
i
z
e
d
W
S
M
S
E
Existing (Orange) and proposed (Blue) algorithms
SNR = 18dB
SNR = 12dB
SNR = 6dB
Fig. 2. Comparison of the proposed algorithm (Algorithm I) and the existing
algorithm for K = 32 and different N.
C. Convergence speed of the proposed algorithm
In this subsection, we study the convergence speed of
the proposed algorithm. In this regard, we set N = 16
and use 3 initialization matrices. The rst is like in Fig.
1, the second is the truncated DFT matrix of size K and
the third initialization is a random matrix taken from the
complex Gaussian distribution with appropriate scaling. For
these initializations, we plot the convergence behavior of the
proposed algorithm in Fig. 3. As can be seen from this gure,
faster convergence speed is achieved when the initialization
is performed as a truncated version of the DFT matrix of
size K. Nevertheless, for all initialization matrices, we have
achieved the same normalized WSMSE in few iterations (i.e.,
n =
40
K
< 2 outer iteration is required for all initializations).
We would like to mention here that we have observed
similar convergence behavior for other settings of g
k
, k,
2
,
N and K.
VI. CONCLUSIONS
In this paper, we propose novel pilot optimization and
channel estimation algorithm for a multiuser massive MIMO
5 10 15 20 25 30 35 40
0.75
0.755
0.76
0.765
0.77
0.775
0.78
0.785
0.79
0.795
0.8
Iteration number
N
o
r
m
a
l
i
z
e
d
W
S
M
S
E
SNR = 0dB
DFT matrix with pilot reuse
Truncated DFT matrix
Random matrix
(a)
5 10 15 20 25 30 35 40
0.67
0.68
0.69
0.7
0.71
0.72
0.73
SNR = 3dB
Iteration number
N
o
r
m
a
l
i
z
e
d
W
S
M
S
E
DFT matrix with pilot reuse
Truncated DFT matrix
Random matrix
(b)
Fig. 3. Convergence speed of the proposed algorithm when N = 16 and
K = 32 [n =
Iteration no
K
40
32
< 2]. (a) SNR = 0dB. (b) SNR = 3dB.
system with K single antenna MSs, arbitrary N pilot symbols
and TDD channel estimation method. The proposed algorithm
is explained as follows. First, we formulate the channel estima-
tion problem as a WSMSE minimization problem containing
pilot symbols and introduced variables. Second, for xed pilot
symbols, the introduced variables are optimized using MMSE
and generalized Rayleigh quotient methods. Finally, for N = 1
and N = K settings, the pilot symbols are optimized using
SDP convex optimization approach, and for the other settings
of N and K, the pilot symbols are optimized by applying
simple iterative algorithm. When N = K, it is shown that
the latter iterative algorithm gives the optimal pilot symbols
achieved by the SDP method. Simulation results conrm that
the proposed algorithm achieves less WSMSE compared to
that of the conventional semi-orthogonal pilot symbol and
MMSE channel estimation algorithm which creates pilot con-
tamination.
REFERENCES
[1] D. P. Palomar, M. A. Lagunas, and J. M. Ciof, Optimum linear joint
transmit-receive processing for MIMO channels with QoS constraints,
IEEE Trans. Sig. Proc., vol. 52, no. 5, pp. 1179 1197, May 2004.
[2] T. E. Bogale, B. K. Chalise, and L. Vandendorpe, Robust transceiver
optimization for downlink multiuser MIMO systems, IEEE Tran. Sig.
Proc., vol. 59, no. 1, pp. 446 453, Jan. 2011.
[3] T. E. Bogale and L. Vandendorpe, Weighted sum rate optimization for
downlink multiuser MIMO coordinated base station systems: Centralized
and distributed algorithms, IEEE Tran. Sig. Proc., vol. 60, no. 4, pp.
1876 1889, Dec. 2011.
[4] D. Tse and P. Viswanath, Fundamentals of Wireless Communication,
Cambridge University Press, 2005.
[5] T. L. Marzetta, Noncooperative cellular wireless with unlimited
numbers of base station antennas, IEEE Tran. Wirel. Comm., vol. 9,
no. 11, pp. 3590 3600, Nov. 2010.
[6] T. L. Marzetta, How much training is required for multiuser mimo?, in
40th Asilomar Conf. on Signals, Systems and Computers, Pacic Grove,
CA, USA, 2006.
[7] Q. N. Hien and E. G. Larsson, EVD-based channel estimation in
multicell multiuser MIMO systems with very large antenna arrays, in
ICASSP, 2012, Kyoto, Japan, 2012, pp. 3249 3252.
[8] S. Boyd and L. Vandenberghe, Convex optimization, Cambridge
University Press, Cambridge, 2004.
[9] T. E. Bogale and L. Vandendorpe, Max-Min SNR signal energy based
spectrum sensing algorithms for cognitive radio networks with noise
variance uncertainty, IEEE Tran. Wirel. comm, 2013.
[10] R. H. Horn and C. R. Johnson, Matrix Analysis, pub-CAMBRIDGE,
CAMBRIDGE, 1985.
[11] D. P. Palomar, J. M. Ciof, and M. A. Lagunas, Joint Tx-Rx beam-
forming design for multicarrier MIMO channels: A unied framework
for convex optimization, IEEE Trans. Sig. Proc., vol. 51, no. 9, pp.
2381 2401, Sept. 2003.