1

Spatial Channel Covariance Estimation
for the Hybrid MIMO Architecture:
A Compressive Sensing Based Approach
arXiv:1711.04207v1 [] 11 Nov 2017

Sungwoo Park and Robert W. Heath Jr.

Abstract
Spatial channel covariance information can replace full knowledge of the entire channel matrix
for designing analog precoders in hybrid multiple-input-multiple-output (MIMO) architecture. Spatial
channel covariance estimation, however, is challenging for the hybrid MIMO architecture because the
estimator operating at baseband can only obtain a lower dimensional pre-combined signal through fewer
radio frequency (RF) chains than antennas. In this paper, we propose two approaches for covariance
estimation based on compressive sensing techniques. One is to apply a time-varying sensing matrix,
and the other is to exploit the prior knowledge that the covariance matrix is Hermitian. We present
the rationale of the two ideas and validate the superiority of the proposed methods by theoretical
analysis and numerical simulations. We conclude the paper by extending the proposed algorithms from
narrowband massive MIMO systems with a single receive antenna to wideband systems with multiple
receive antennas.

I. I NTRODUCTION
One way to increase coverage and capacity in wireless communication systems is to use a
large number of antennas. For instance, millimeter wave systems use large antenna arrays to
obtain high array gain, thereby increasing cellular coverage [1], [2]. Sub-6 GHz systems are
also likely to equip many antennas at a base station (BS) to increase cellular spectral efficiency
by transmitting data to many users simultaneously via massive MIMO systems [3], [4]. Using
a large number of antennas in a conventional MIMO architecture, however, results in high cost
S. Park and R. W. Heath Jr. are with the Wireless Networking and Communication Group (WNCG), Department of Electrical
and Computer Engineering, The University of Texas at Austin, TX, 78701 USA. (e-mail: {swpark96, rheath}@utexas.edu).
This work is supported in part by the National Science Foundation under Grant No. 1711702, and by a gift from Huawei
Technologies.

2

and high power consumption because each antenna requires its own RF chain [5]. To solve
this problem, hybrid analog/digital precoding reduces the number of RF chains by dividing the
linear process between the analog RF part and the digital baseband part. Several hybrid precoding
techniques have been proposed for single-user MIMO [6]–[9] and multi-user MIMO [10]–[14]
in either sub-6 GHz systems or millimeter wave systems.
Most prior work on hybrid precoding assumes that full channel state information at the
transmitter (CSIT) for all antennas is available when designing the analog precoder. Full CSIT,
however, is difficult to obtain even in time-division duplexing (TDD) systems if many antennas are employed. Furthermore, the hybrid structure makes it more difficult to estimate the
full CSIT of the entire channel matrix for all antennas because the estimator can only see
a lower dimensional representation of the entire channel due to the reduced number of RF
chains. Although some channel estimation techniques for the hybrid architecture were proposed
by using compressive sensing techniques that exploit spatial channel sparsity [15]–[17], these
techniques assume that channel does not change during the estimation process which requires
many measurements over time. Consequently, these channel estimation techniques can be applied
only to slowly varying channels.
Unlike the hybrid precoding techniques that require full CSIT, there exists another type of
hybrid precoding that uses long-term channel statistics such as spatial channel covariance instead
of full CSIT in the analog precoder design [13], [14], [18], [19]. This approach has advantages
over using full CSIT. First, the spatial channel covariance is well modeled as constant over
many channel coherence intervals [18], [20]. Second, the spatial channel covariance is constant
over frequency in general [21], so covariance information is suitable for the wideband hybrid
precoder design problem where one common analog precoder must be shared by all subcarriers
in wideband OFDM systems [9]. For these reasons, spatial channel covariance information is a
promising alternative to full CSIT for hybrid MIMO architecture.
Although several techniques have been developed to estimate spatial channel covariance or
its simplified version such as subspace or angle-of-arrival (AoA) for hybrid architecture, these
methods have practical issues. For example, the sparse array developed in [22] and the coprime
sampling used in [23] reduce the number of RF chains for the covariance or AoA estimation
as they disregard some of the antennas by exploiting redundancy in linear arrays of equidistant
antennas. These methods, however, have a limitation on the configuration of the number of
RF chains and antennas. In [24], various estimation methods were proposed based on convex

3

optimization problems coupled with the coprime sampling. The proposed methods, though,
require many iterations to converge. A subspace estimation method using the Arnoldi iteration
was proposed for millimeter wave systems in [25]. The estimation target of this method is not
the subspace of a spatial channel covariance matrix but that of a channel matrix itself, and
the channel matrix must be constant over time during the estimation process. Consequently,
it is difficult to apply the method in [25] to time-varying channels. The subspace estimation
method proposed in [26] does not use any symmetry properties of the covariance matrix, while
the AoA estimation for the hybrid structure in [27] does not exploit the sparse property of
millimeter wave channels. In [28], [29], AoA estimation methods for the hybrid structure were
proposed by applying compressive sensing techniques with vectorization but do not fully exploit
the slow varying AoA property. Instead of using a vector-type compressive sensing, matrixtype compressive sensing techniques were developed for so-called multiple measurement vector
(MMV) problems [30]–[33]. Most prior work on the MMV problems, however, assumes that
the sensing matrix is fixed over time to model the problem in a matrix form, which is not an
efficient strategy when the measurement size becomes as small as the sparsity level.
In this paper, we propose a novel spatial channel covariance estimation technique for the
hybrid MIMO architecture. Considering spatially sparse channels, we develop the estimation
technique based on compressive sensing techniques that leverage the channel sparsity. Between
two different approaches in the compressive sensing field (convex optimization algorithms vs.
greedy algorithms), we focus on the greedy approach because they provide a similar performance
to convex optimization algorithms in spite of their simpler complexity [34]. Based on well-known
greedy algorithms such as orthogonal matching pursuit (OMP) and simultaneous OMP (SOMP),
we improve performance by applying two key ideas: one is to use a time-varying sensing matrix,
and the other is to exploit the fact that covariance matrices are Hermitian. Motivated by the fact
that SOMP is a generalized version of OMP, we develop a further generalized version of SOMP
by incorporating the proposed ideas. The algorithm modification is simple, but the performance
improvement is significant, which is demonstrated by numerical and analytical results.
We first develop the spatial channel covariance estimation work for a simple scenario where
a mobile station (MS) has a single antenna in narrowband systems. Preliminary results were
presented in the conference version of this paper [35]. In this paper, we add two new contributions
to our prior work. First, we present theoretical analysis to validate the rationale of the two
proposed ideas and show the superiority of the proposed algorithm. The theoretical analysis

4

demonstrates that the use of a time-varying sensing matrix dramatically improves the covariance
estimation performance, in particular, when the number of RF chains is not so large and similar
to the number of channel paths. The analytical results also disclose that exploiting the Hermitian
property of the covariance matrix provides an additional gain. Second, we extend the estimation
work to other scenarios. Considering that an MS as well as a BS has hybrid architecture with
multiple antennas and RF chains, we modify the proposed algorithm to adapt to the situation
where analog precoders as well as analog combiners change over time. We also modify the
algorithm for the wideband systems by using the fact that frequency selective baseband combiners
do not improve the estimation performance.
We use the following notation throughout this paper: A is a matrix, a is a vector, a is a
scalar, and A is a set. kak0 and kak2 are the l0 -norm and l2 -norm of a vector a. kAkF denotes
a Frobenius norm. AT , AC , A∗ , and A† are transpose, conjugate, conjugate transpose, and
Moore-Penrose pseudoinverse. [A]i,: and [A]:,j are the i-th row and j-th column of the matrix
A. [A]:,S denotes a matrix whose columns are composed of [A]:,j for j ∈ S. If there is nothing
ambiguous in the context, [A]:,j , [A]:,S , [At ]:,j , and [At ]:,S are replaced by aj , AS , at,j , and
At,S . A \ B is the relative complement of A in B, diag(A) is a column vector whose elements
are the diagonal elements of A, and |A| is the cardinality of A. A ⊗ B, A } B, and A

B

denote the Kronecker product, the Hadamard product, and the column-wise Khatri-Rao product.
(K)

A
(K)

A

B is the generalized Khatri-Rao product with respect to K partitions, which is defined as
i
h
i
i
h
h
B = A1 ⊗ B1 · · · AK ⊗ BK where A = A1 · · · AK and B = B1 · · · BK .

II. S YSTEM MODEL AND PRELIMINAIRES : COVARIANCE ESTIMATION VIA COMPRESSIVE
SENSING BASED CHANNEL ESTIMATION

In this section, we present the system model and briefly overview prior work on compressive
sensing based channel estimation followed by covariance calculation.
A. System model
Consider a system where a BS with N antennas and M RF chains communicates with an
MS. We focus on a narrowband massive MIMO system where an MS with a single antenna; we
extend the work to the multiple-MS-antenna case and the wideband case in Section V. Let L be
the number of channel paths, g`,t be a time-varying channel coefficient of the `-th channel path
at the t-th snapshot, and a(φ` ) be an array response vector associated with the AoA of the `-th

5

channel path φ` . Assuming that AoAs do not change during T snapshots, the uplink channel at
the t-th snapshot can be represented as
ht =

L
X

g`,t a(φ` ).

(1)

`=1

Considering spatially sparse channels, the channel model in (1) can be approximated in the
compressive sensing framework as
ht = Agt ,

(2)

where A ∈ CN ×D is an dictionary matrix whose D columns are composed of the array response
vectors associated with a predefined set of AoAs, and gt ∈ CD×1 is a sparse vector with only
L nonzero elements whose positions and values correspond to their AoAs and path gains.
Let dt be a training symbol with |dt | = 1, and zt be a Gaussian noise with CN (0, σ 2 I). The
received signal can be represented as
xt = Agt dt + zt .

(3)

Combined with an analog combining matrix Wt ∈ CM ×N , the received signal multiplied by d∗t
at baseband is expressed as
yt = d∗t Wt xt = Φt gt + Wt nt ,

(4)

where Φt = Wt A denotes an overall sensing matrix and nt = d∗t zt ∼ CN (0, σ 2 I). By letting
Rg = E[gt gt∗ ], the spatial channel covariance matrix is expressed as Rh = E[ht h∗t ] = ARg A∗ .
The goal is to estimate Rh with given y1 , ..., yT .
B. Preliminaires: covariance estimation via compressive sensing based channel estimation
Some prior work developed compressive channel estimators for the hybrid architecture [15]–
[17]. Once the channel vectors are estimated, the covariance can be calculated from the estimated
channels. In this section, we overview two estimation approaches to make a comparison.
A baseline technique to estimate the channel vector at each time is formulated as
min kyt − Φt gt k2 s.t. kgt k0 ≤ L.
gt

(5)

This optimization problem is known as a single measurement vector (SMV) problem and can

6

Algorithm 1 OMP

Algorithm 2 SOMP

Input: Φ, y, and L
Initialize: v = y, S = Ø, ĝ = 0
for n = 1 : L do
j = arg maxi |φ∗i v|
S =
S ∪ {j}

v = I − ΦS Φ†S y
end for
ĝS = Φ†S y
Output: ĝ

Input: Φ, Y, and L
Initialize: V = Y, S = Ø, Ĝ = 0
for n = 1 : L do
j = arg maxi kφ∗i Vk2
S = S ∪ {j}

V = I − ΦS Φ†S Y
end for
[Ĝ]S,: = Φ†S Y
Output: Ĝ

be solved by the OMP algorithm described in Algorithm 1, which is a simple but efficient
approach among various solutions to the SMV problem [36]. Once the gt ’s are estimated during
T snapshots, then the covariance matrix of the channel can be calculated as
!
T
1X
R̂h = A
gt gt∗ A∗ .
T t=1

(6)

In contrast to SMV, another compressive sensing technique called MMV [30] exploits the
observation that g1 , ..., gT shares a common support if the AoAs do not change during the
estimation process. The received signal model in (4) can be written in a matrix form as
Y = ΦG + WN,

(7)

h
i
h
i
h
i
where Y = y1 · · · yT , G = g1 · · · gT , and N = n1 · · · nT . Note that in this
MMV problem format, the sensing matrix Φ must be common over time. The optimization
problem in the MMV is formulated as
min kY − ΦGkF s.t. kGkrow,0 ≤ L,
G

(8)

S
where kGkrow,0 represents the row sparsity defined as kGkrow,0 = | t supp ([G]:,t )|. This
optimization problem can be solved by the SOMP algorithm described in Algorithm 2 [37],
[38]. Regarding the selection rule of SOMP, the `2 -norm in Algorithm 2 can be replaced by the
`1 -norm because there is no significant difference between the two in terms of performance [31].
We use `2 -norm throughout this paper for analytical tractability. Although SOMP is known to
outperform OMP in general, SOMP still needs improvement, particularly when the number of

7

RF chains is small, as will be discussed in Section III-A.
III. C OVARIANCE ESTIMATION FOR HYBRID ARCHITECTURE : TWO KEY IDEAS
In this section, we present two key ideas to improve covariance estimation work based
on compressive sensing techniques. First, we develop an estimation algorithm by applying a
time-varying sensing matrix. Second, we propose another algorithm that exploits the Hermitian
property of the covariance matrix. Finally, we combine the two proposed algorithms.
A. Applying a time-varying analog combining matrix
Although compressive sensing can reduce the required number of measurements, there is a lim
itation on reducing the measurement size. This lower bound is known to be M = O L log D
L
for both SMV and MMV when the sensing matrix satisfies the restricted isometry property (RIP)
condition. [39]. One possible solution to overcome this limitation is to extend the measurement
vector size into the time domain, i.e., gathering the measurement vectors of multiple snapshots
with different sensing matrices over time. The same approach is also used in [15] and [36],
where the measurements are stacked together as

  

  
Φ
W1 n1
W1 A
y1

  1

  
 .   . 
 ..   .. 
 .  =  .  g +  ..  =  ..  g + ñ.

  

  
ΦT
WT nT
WT A
yT

(9)

Note that the key assumption in (9) is that g is constant during the estimation process. Consequently, this technique can be applied only to static channels, not time-varying channels.
It is worthwhile to compare the MMV signal model in (7) with the signal model where the
time-varying combining matrix is used for the static channel as shown in (9). In (7), yt ’s are
stacked in columns due to the fact that gt changes over time but Φ is fixed. In contrast, the
signal model in (9) adopts a row-wise stack of yt because Φt changes over time but g is fixed.
If both Φt and gt change over time, we can not stack yt ’s in either columns or rows and thus
need another approach. Regarding this scenario, two questions arise: 1) is it useful to employ a
time-varying sensing matrix for a time-varying channel? 2) how can we recover the data in this
time-varying sensing matrix and time-varying data? We will show that the time-varying sensing
matrix can increase the recovery success rate especially when M is not much larger than L by

8

analysis in Section IV and simulations in Section VI. In this subsection, we develop a recovery
algorithm to answer the problem of recovery with a time-varying sensing matrix.
To apply the time-varying matrix concept to SOMP, we first focus on the fact that SOMP is
a generalized version of OMP. The reason is that the optimization problem of MMV in (8) can
be rewritten as
min

g1 ,...,gT

T
X

kyt − Φgt k22 s.t.

T
[

supp (gt ) ≤ L.

(10)

t=1

t=1

Note that SOMP is equivalent to OMP when T = 1. Unlike the original formulation in (8), this
reformulated form gives an insight into how to apply the time-varying sensing matrix. We can
further generalize the optimization problem in (10) by replacing Φ with Φt to reflect the timevarying feature of the sensing matrix. Noting that the selection rule of SOMP in Algorithm 2 is
P
equivalent to j = arg maxi kφ∗i Vk22 = arg maxi Tt=1 |φ∗i vt |2 , we modify the selection rule by
replacing φi with φt,i (= [Φt ]:,i ) to adapt to the time-varying sensing matrix case. The residual




matrix in SOMP, V = I − ΦS Φ†S Y, can also be replaced by vt = I − Φt,S Φ†t,S yt for
t = 1, ..., T . The modified version of SOMP, which we call dynamic SOMP (DSOMP), is
described in Algorithm 3.
B. Exploiting the Hermitian property of the covariance matrix
The covariance estimation techniques using OMP, SOMP, and DSOMP in the previous subsections employ a two-step approach where the channel gain vectors gt ’s (and thus channel
vectors ht ’s) are estimated and then the covariance matrix is calculated from (6). If it is not
the channel but the covariance that needs to be estimated, the first step estimating the channel
explicitly is unnecessary. In this subsection, we take a different approach that directly estimates
the covariance Rg without estimating the instantaneous channel gain vectors.
The relationship between Rg and Ry is given by
Ry = ΦRg Φ∗ + σ 2 WW∗ .

(11)

If we consider a special case where Rg is assumed to be a sparse diagonal matrix as in Fig. 1(c),
the relationship between Rg and Ry in (11) can be rewritten as
vec(Ry ) = ΦC


Φ diag(Rg ) + σ 2 vec(WW∗ ),

(12)

9

(a) Unstructured

(b) Sparse rows

(c) Sparse diagonal (d) Sparse Hermitian

Fig. 1. Sparse matrix types: (a) an unstructured sparse matrix where the elements are randomly spread, (b) a structured sparse
matrix with sparse rows, (c) a structured sparse diagonal matrix, and (d) a structured sparse Hermitian matrix.

by using vectorization and the column-wise Khatri-Rao product, then OMP can be directly
applied to this reformulated problem without any modification [40]. This approach, however,
has a limitation in application to realistic scenarios because the covariance Rg is not a diagonal
matrix in general [35]. Instead of assuming that Rg is a sparse diagonal matrix, we consider
Rg to be a sparse Hermitian matrix as in Fig. 1(d). As MMV uses the fact that G is a matrix
with sparse rows as in Fig. 1(b), we exploits the Hermitian property of the covariance matrix
Rg . The optimization problem is represented in a matrix form like MMV,
min kRY − ΦRg Φ∗ kF s.t. kRg klattice,0 ≤ L,

(13)

Rg

where kRg klattice,0 is defined as kRg klattice,0 =

S

i

supp ([Rg ]:,i )

S

j

supp ([Rg ]j,: ) .

A greedy approach similar to OMP and SOMP can be applied to find a suboptimal solution.
At each iteration, the algorithm needs to find the solution to the following sub-problems as
min kRY − ΦSn RX Φ∗Sn kF ,

(14)

RX

where RX is a Hermitian matrix whose dimension is less than or equal to L. By using the least
squares method along with vectorization, the optimal solution to (14) is given by
(opt)

RX
where Φ†Sn = Φ∗Sn ΦSn

−1


∗
= Φ†Sn RY Φ†Sn ,

(15)
(opt)

Φ∗Sn . Note that if RY is Hermitian, RX

is also Hermitian. The

proposed algorithm using (15), which we call covariance OMP (COMP), is described in Algorithm 4. Note that the proposed COMP uses quadratic forms instead of linear forms.
The time-varying sensing matrix concept in Section III-A can also be applied to COMP. The
generalization of COMP, however, requires a careful modification. In covariance estimation based

10

Algorithm 3 Dynamic SOMP (DSOMP)

Algorithm 4 Covariance OMP (COMP)

Input: Φ1 , ..., ΦT , Y, L
Initialize: V = Y, S = Ø, Ĝ = 0
for n = 1 : L doP
2
j = arg maxi Tt=1 φ∗t,i vt
S = S ∪ {j}

†
vt = I − Φt,S Φt,S yt , ∀t
end for h
i
[Ĝ]S,: = Φ†1,S y1 · · · Φ†T,S yT

Input: Φ, Y, L
Initialize: V = R̂y = T1 YY∗ , S = Ø, R̂g = 0
for n = 1 : L do
j = arg maxi φ∗i Vφi
S = S ∪ {j}
∗

†
†
V = R̂y − ΦS ΦS R̂y ΦS ΦS
end for
 ∗
[R̂g ]S,S = Φ†S R̂y Φ†S

Output: Ĝ

Output: R̂g

Algorithm 5 Dynamic COMP (DCOMP)
Input: Φ1 , ..., ΦT , Y, L
Initialize: Vt = R̂y,t = yt yt∗ , ∀t, S = Ø, R̂g = 0
for n = 1 : L doPT
j = arg maxi t=1 φ∗t,i Vt φt,i
S = S ∪ {j}
∗

Vt = R̂y,t − Φt,S Φ†t,S R̂y,t Φt,S Φ†t,S , ∀t
end for

∗
PT
†
†
1
[R̂g ]S,S = T t=1 Φt,S R̂y,t Φt,S
Output: R̂g

on COMP, the calculation of R̂y is followed by the estimation of R̂g . This is because R̂y can
be represented as a function of R̂g as
T
1X
yt yt∗ = Φ
R̂y =
T t=1

!
T
1X
gt gt∗ Φ∗ = ΦR̂g Φ∗ .
T t=1

(16)

If the sensing matrix changes over time, R̂y can not be represented as a function of R̂g because
T
T
1X
1X
∗
R̂y =
y t yt =
Φt gt gt∗ Φ∗t .
T t=1
T t=1

(17)

For this reason, instead of using the sample covariance matrix of the measurement vectors
R̂y , we use a per-snapshot yt yt∗ , which can be regarded as a one sample estimate of the sample
covariance matrix. Note that this extreme sample covariance is not a diagonal matrix in general
even when the channel paths are uncorrelated. Using the fact that yt yt∗ for all snapshots are sparse
Hermitian matrices sharing the same positions of nonzero elements, we develop the dynamic

11

COMP (DCOMP) from COMP in a similar way that we develop DSOMP from SOMP, which is
described in Algorithm 5. Note that DCOMP becomes equivalent to COMP if Φ1 = · · · = ΦT .
IV. T HEORETICAL ANALYSIS OF THE PROPOSED ALOGIRHTMS
In this section, we analyze the benefit of the time-varying sensing matrix and compare SOMP,
DSOMP, and DCOMP. In summary, the recovery success probability of SOMP, DSOMP, and
DCOMP increases as the number of measurements T increases. The recovery success probability
of SOMP, however, saturates and does not approach one as T goes to infinity even in the noiseless
case if the number of RF chains M is not much larger than the number of channel paths L.
In contrast to SOMP, the recovery success probability of DSOMP and DCOMP approaches one
as T increases for any M and L. In other words, DSOMP and DCOMP can guarantee perfect
recovery even when the number of RF chains is so small that both OMP and SOMP fail to
recover the support even in the noiseless case. Finally, we will show that DCOMP has a higher
recovery success probability than DSOMP.
Let us first consider how to design the sensing matrix Φ. For analytical tractability, we confine
the dictionary size to D = N in this section. The algorithms, however, can be applied to more
h
i
general cases such as D > N . The sensing matrix can be represented as Φ = φ1 · · · φN =
WA. One possible option for the sensing matrix is to use a random analog precoding matrix
W = WRF such that each element in WRF has a constant amplitude and a random phase with
an independent uniform distribution in [0, 2π]. Since the elements in the sensing matrix designed
in this way have a sub-Gaussian distribution, this simple random sensing matrix Φ satisfies the
RIP condition [34]. There are other desirable features that the sensing matrix needs to have.
For example, it is desirable for the sensing matrix to have a small mutual coherence defined as
q
|φ∗j φi |
N −M
[34].
The
mutual
coherence
has
a
lower
bound
, which is known
ρ = maxi6=j kφ kkφ
k
M (N −1)
j

i

as Welch bound. This bound can be achievable if the column vectors in Φ constitute an equalnorm equiangular tight frame, i.e., Φ satisfies that 1) kφi k = c1 , ∀i, 2) |φ∗j φi | = c2 , ∀j 6= i, and
3) ΦΦ∗ = c3 I. It is impossible for a sensing matrix Φ that is designed with a random phase
analog combining matrix WRF to meet these three conditions at all times. Nevertheless, it is
possible to meet the tight frame condition for any random analog precoding matrix WRF if the
baseband combining matrix WBB is employed such that
1

∗
W = WBB WRF = (WRF WRF
)− 2 WRF .

(18)

12

If the baseband combining matrix is designed as (18), the columns in the overall sensing matrix
Φ constitute the tight frame because it satisfies
∗
∗
ΦΦ∗ = WBB WRF AA∗ WRF
WBB
= N I.

(19)

Consequently, any sensing matrix can be transformed to satisfy the tight frame condition by
1

∗
using WBB = (WRF WRF
)− 2 .

The following theorem provides the exact condition for the perfect recovery condition at each
iteration of the DSOMP algorithm. Note that OMP and SOMP are special cases of DSOMP.
Theorem 1 Let S ⊂ N = {1, 2, ..., N } be an optimal support set with |S| = L and G be an
ideal channel gain matrix. Suppose that So ⊂ S with |So | = Lo was perfectly recovered at the
previous Lo iterations in the DSOMP algorithm. In the noiseless case, one of the elements in
the optimal support S can be perfectly recovered at the (Lo + 1)-th iteration if and only if
ρ(DS) (S, So , Φ1 , ..., ΦT , G) =



where Ψt = I −

Φt,So Φ†t,So



maxj∈N \S

PT

P

maxj∈S\So

PT

P

t=1
t=1

i∈S\So

ψ ∗t,j ψ t,i gt,i

∗
i∈S\So ψ t,j ψ t,i gt,i

2
2

< 1,

(20)

Φt .

Proof: See Appendix A.
Although the newly defined metric ρ(DS) (S, So , Φ1 , ..., ΦT , G) identifies the exact condition
for the perfect recovery condition, its dependence on the input data G makes this metric less
attractive. Some prior work takes an approach to find upper bounds of certain coherence-related
metrics for any random input data [41], but the bounds are loose, which can be regarded as a
worst case analysis. In addition, those upper bounds are meaningful only for the case where M
is significantly larger than L. In the extreme case where M and L are similar in value, the upper
bounds are too loose to give a useful insight.
Instead of an upper bound, we focus on the distribution of ρ(DS) (S, So , Φ1 , ..., ΦT , G) considering S, So , Φ1 , ..., ΦT , and G as random variables. While the exact distribution depends on the
distribution of G, Proposition 1 and Theorem 2 show that ρ(DS) (S, So , Φ1 , ..., ΦT , G) converges
to a value that does not depend on G as T becomes larger.

13

Proposition 1 Suppose that gt,i are independent random variables with zero mean and unit
variance. Let ρ(S) (S, So , Φ, G) denote ρ(DS) (S, So , Φ1 , ..., ΦT , G) in the SOMP case where Φ1 =
· · · = ΦT = Φ. As T → ∞, ρ(S) (S, So , Φ, G) converges to
(S)

ρ (S, So , Φ, G) →

(S)
ρT →∞ (S, So , Φ)

=

maxj∈N \S

P

ψ ∗j ψ i

maxj∈S\So

P

ψ ∗j ψ i

i∈S\So
i∈S\So

2
2,

(21)



where Ψ = I − ΦSo Φ†So Φ.
Proof: See Appendix B.
(S)

Note that ρT →∞ (S, So , Φ) in (21) can be rewritten as

P
2
∗
maxj∈N \S
i∈S\So ψ j ψ i
(S)

.
ρT →∞ (S, So , Φ) =
P
2
∗
maxj∈S\So |ψ j |4 +
ψj ψi
i6=j

(22)

i∈S\So

In the OMP case, ρ(DS) (S, So , Φ1 , ..., ΦT , G) becomes equivalent to
P

2
∗
maxj∈N \S
ψ
ψ
+
X
j
j i
i∈S\So

,
ρ(O) (S, So , Φ, g) =
P
2
∗
ψ j ψ i + Xj
maxj∈S\So |ψ j |4 +
i6=j

(23)

i∈S\So

where Xj =

P

i1 ∈S\So

P

i2 ∈S\So
i2 6=i1

gi1 gi∗2 ψ ∗j ψ i1 ψ ∗i2 ψ j ∈ R. Since the addition of a random variable

Xj results in increasing the variance of the metric at j, the maximum values in both the numerator
and the denominator for the OMP case are likely to be larger than those of the SOMP case.
Compared with the SOMP case, the relative rate of increase in the numerator is higher than that
of increase in the denominator with high possibility because |ψ j |4 is larger than |ψ ∗j ψ i |2 for
i 6= j in general for a random sensing matrix. Consequently, we can infer that the normalized
coherence metric of SOMP is likely to be less than that of OMP. We validate this inference by
simulations for random Φ, S, and So in Section VI-A.
Although this explains the superiority of SOMP over OMP, we can see that the probability of
(S)

ρT →∞ (S, So , Φ) < 1 does not become one in general. This means that SOMP can not guarantee
the perfect recovery even when an infinite number of snapshots are measured for the recovery.
Unlike SOMP, DSOMP can always guarantee perfect recovery as T becomes larger as shown
in the following theorem.

14

Theorem 2 As T → ∞, ρ(DS) (S, So , Φ1 , ..., ΦT , G) in (20) converges to a deterministic value
(DS)

ρT →∞ (N, M, L, Lo ) that has an upper bound as
(DS)

ρ(DS) (S, So , Φ1 , ..., ΦT , G) → ρT →∞ (N, M, L, Lo )
−1

M − N + (M − Lo )(N − Lo − 1)
,
≤ 1+
(N − M )(L − Lo )

(24)

and this upper bound is always less than or equal to one for any N ≥ M ≥ L > Lo , i.e,. perfect
recovery is guaranteed.
Proof: See Appendix C.
In Section VI, we show by simulations that a reasonably small T can provide a significant
gain over SOMP even in the extreme case when M = L. In addition, the converged value
(DS)

of DSOMP ρT →∞ (N, M, L, Lo ) is not a function of sensing matrices Φ1 , ..., ΦT while that of
(S)

SOMP ρT →∞ (S, So , Φ) depends on the sensing matrix Φ. In other words, in contrast to the
SOMP case where sophisticated sensing matrix design is crucial to the estimation performance,
DSOMP can reduce the necessity of the elaborate design of the sensing matrix as T increases.
Similar to the DSOMP case, DCOMP can also guarantee the perfect recovery. The perfect
recovery condition of the DCOMP case is shown in the following theorem.
Theorem 3 Let S be an optimal support set with |S| = L and G be an ideal channel gain
matrix. Suppose that So ⊂ S with |So | = Lo was perfectly recovered at the previous Lo iterations
in the DCOMP algorithm. In the noiseless case, one of the elements in the optimal support S
can be perfectly recovered at the (Lo + 1)-th iteration if and only if
ρ

(DC)

(S, So , Φ1 , ..., ΦT , G) =

(DC)
∗
t=1 φt,j Qt,S,So φt,j
P
(DC)
maxj∈S\So Tt=1 φ∗t,j Qt,S,So φt,j

maxj∈N \S

PT

< 1,

(25)

∗
∗
where Qt,S,So = Φt,S gt,S gt,S
Φ∗t,S − Pt,So Φt,S gt,S gt,S
Φ∗t,S Pt,So and Pt,So = Φt,So Φ†t,So .
(DC)

Proof: From the definition of the residual matrix Vt in Algorithm 5, the selection rule of
DCOMP can be represented as
j

(opt)

= arg max

and this completes the proof.

j∈N

T
X
t=1

φ∗t,j Vt φt,j

= arg max
j∈N

T
X
t=1

(DC)

φ∗t,j Qt,S,So φt,j ,

(26)

15

(DS)

(DC)

∗
Let Qt,S,So = (I − Pt,So ) Φt,S gt,S gt,S
Φ∗t,S (I − Pt,So ). Notice that, if Qt,S,So is replaced by
(DS)

Qt,S,So , then Theorem 3 becomes identical to Theorem 1. It is also worthwhile to note that both
(DS)

(DC)

DSOMP and DCOMP select the same support at the first iteration because Qt,S,So = Qt,S,So for
So = ∅. The superiority of DCOMP over DSOMP starts from the second iteration, which can
be inferred from the following theorem.
Theorem 4 ρ(DC) (S, So , Φ1 , ..., ΦT , G) in (25) converges to a deterministic value as T → ∞,
(DC)

and the deterministic value ρT →∞ (N, M, L, Lo ) always satisfies
(DC)

(DS)

ρT →∞ (N, M, L, Lo ) ≤ ρT →∞ (N, M, L, Lo ),

(27)

for any N ≥ M ≥ L > Lo .
Proof: See Appendix D.
Until now, we considered a time-varying analog combining matrix and its associated digital
combining matrix. One question arises: if the analog combiner is fixed and only the digital
combiner is time-varying, does the time variance of the baseband combiner improve the recovery
performance? The following proposition shows that the time-varying baseband combining matrix
itself does not improve the performance when the analog combining matrix is fixed.
1

∗
)− 2 for a fixed WRF and a random unitary maProposition 2 Let WBB,t = Ut (WRF WRF

trix Ut ∈ CM ×M such that Φt = WBB,t WRF A constitutes a tight frame for any t. Then,
ρ(DS) (S, So , Φ1 , ..., ΦT , G) of DSOMP becomes identical to ρ(DS) (S, So , Φ0 , ..., Φ0 , G) of SOMP
1

∗
)− 2 WRF A.
where Φ0 = (WRF WRF

Proof: Each term in both the numerator and the denominator in (51) becomes
(DS)

∗
φ∗t,j Qt,S,So φt,j = φ∗t,j (I − Pt,So ) Φt,S gt,S gt,S
Φ∗t,S (I − Pt,So ) φt,j

(28)

∗
= φ∗0,j (I − P0,So ) Φ0,S gt,S gt,S
Φ∗0,S (I − P0,So ) φ0,j ,

because Φt = Ut Φ0 and Pt,So = Ut P0,So U∗t for all t. Since this is the same as in the SOMP
case where Φ1 = · · · = ΦT = Φ0 , the proof is completed.
The result in Proposition 2 is beneficial especially to the wideband case where the analog
combining matrix must be fixed over frequency. According to Proposition 2, we will use the
frequency-invariant baseband combining matrix when modifying the algorithm for the wideband
case. This will be discussed in more detail in Section V-B.

16

frame t

frame t-1

frame t+1

MT training symbols data
(s = 1, 2, …, MT)
Fig. 2. Frame structure in the case when an MS has hybrid structure with multiple RF chains and antenna. Consecutive MT
training symbols are inserted in every frame, and each training symbol is transmitted one at a time via its dedicated RF chain.

V. E XTENSION TO OTHER SCENARIOS
In Section III, we considered a narrowband system and assumed that an MS has a single
antenna. In this section, we first extend the work to the multiple-MS-antenna case and then to
the wideband case.
A. Extension to the multiple-MS-antenna case
In this subsection, we extend the work in Section III to the case where the MS has hybrid
architecture with multiple antennas and RF chains. The signal model is shown in Fig. 2 where
NT , MT , NR , and MR denote the number of transmit antennas, transmit RF chains, receive
antennas, and receive RF chains. We assume that consecutive MT training symbols are used in
each frame as shown in Fig. 2, and the symbol index per frame is denoted by s. We assume
that the duration of MT symbols is less than the channel coherence time, i.e., the channel is
invariant during MT consecutive symbol transmission. Let θ` be the angle of departure (AoD)
of the `-th path. The channel model for the single-MS-antenna case in (1) can be extended to
the multiple-MS-antenna case as
Ht =

L
X

gt,` aR (φ` )a∗T (θ` ).

(29)

`=1

Similarly to (2), this can also be represented in the compressive sensing framework as Ht =
AR Gt A∗T , where AR ∈ CNR ×DR and AT ∈ CNT ×DT are dictionary matrices associated with
AoA and AoA, and Gt ∈ CDR ×DT is a sparse matrix with L nonzero elements with associated
with complex channel path gains.
At each symbol s at frame t, one training signal is transmitted through one transmit RF chain.
Let Wt,s be a combiner and ft,s a precoder at frame t and symbol s. Assuming that the training

17

symbol is known to the BS, thereby omitted as in the single-MS-antenna case, the received
signal at baseband is represented as
yt,s = Wt,s AR Gt A∗T ft,s + Wt,s nt,s .

(30)

There are four different types with respect to how to apply Wt,s and ft,s within a frame: 1)
fixed Wt,s and fixed ft,s , 2) time-varying Wt,s and fixed ft,s , 3) fixed Wt,s and time-varying
ft,s , and 4) time-varying Wt,s and time-varying ft,s .
In the first approach, all the precoders are the same within a frame, and so are the combiners
such that ft,1 = · · · = ft,MT = ft and Wt,1 = · · · = Wt,MT = Wt . Since these combiners
and precoders are constant within a frame in this approach, averaging the received signals at
baseband within a frame can increase the effective signal-to-noise ratio (SNR). The averaged
received signal is represented as
MT
1 X
yt,s = Wt AR Gt A∗T ft + Wt n̄t ,
ȳt =
MT s=1

where n̄t =

1
MT

(31)



σ2
n
.
Note
that
n̄
∼
CN
0,
I
, which means that this process
t
s=1 t,s
MT

PM T

averages out the noise, increasing the effective SNR. The signal model per frame in (31) can be
transformed into
ȳt = ftT ⊗ Wt




AC
T ⊗ AR vec(Gt ) + Wt n̄t .

(32)

In the second approach, where only combiners vary and precoders are fixed within a frame
such that ft,1 = · · · = ft,MT = ft , the first step is to stack yt,1 , ..., yt,MT in rows. This row-wise
stack yields a MR MT × 1 vector per frame,
ỹt,agg = W̃t,agg AR Gt A∗T ft + ñt,agg ,
where W̃t,agg

(33)

h
iT
h
iT
T
T
T
T
T
T
= Wt,1 · · · Wt,MT and ñt,agg = nt,1 Wt,1 · · · nt,MT Wt,MT . The signal

model in (33) can be rewritten as



ỹt,agg = ftT ⊗ W̃t,agg AC
T ⊗ AR vec(Gt ) + ñt,agg .

(34)

The third approach changes only precoders, and combiners are fixed within a frame such that

18

Algorithm 6 Wideband DCOMP (WB-DCOMP)
Input: Φ1 , ..., ΦT , y1,1 , ..., y
T,K , L
P
∗
Initialize: Vt = R̂y,t = K1 K
k=1 yt,k yt,k , S = Ø, R̂g = 0
for n = 1 : L doP
j = arg maxi Tt=1 φ∗t,i Vt φt,i
S = S ∪ {j}

∗
Vt = Ry,t − Φt,S Φ†t,S Ry,t Φt,S Φ†t,S , ∀t
end for

∗
PT
†
†
1
[R̂g ]S,S = T t=1 Φt,S Ry,t Φt,S
Output: R̂g

Wt,1 = · · · = Wt,MT = Wt . In contrast to the second approach, yt,1 , ..., yt,MT are stacked in
columns. This column-wise stack makes the aggregate signal model per frame as
Ỹt,agg = Wt AR GD,t A∗T F̃t,agg + Ñt,agg ,

(35)

h
i
h
i
h
i
where Ỹt,agg = yt,1 · · · yt,MT , F̃t,agg = ft,1 · · · ft,MT , and Ñt,agg = Wt∗ nt,1 · · · nt,MT .
By using vectorization, the matrix Ỹt,agg in (35) can be transformed into a vector form as



ỹt,agg = vec(Ỹt,agg ) = F̃Tt,agg ⊗ Wt AC
T ⊗ AR vec(Gt ) + vec(Ñt,agg ).

(36)

For the last approach, both precoders and combiners change within a frame. Like the second
approach, we stack yt,1 , ..., yt,MT in rows. Then, the row-wise stacked vector ỹt,agg becomes

ỹt,agg =

F̃t,agg

(MT )

T
W̃t,agg

T


AC
T ⊗ AR vec(Gt ) + ñt,agg .

(37)

Note that the signal models for all four different approaches in (32), (34), (36), and (37) can be
generalized as ỹt,agg = Θt,agg Aagg gt,agg + ñtt,agg where gt,agg = vec(Gt ), Aagg = AC
T ⊗ AR , and
Θt,agg are dependent on each approach. This generalized signal model has the same format as that
of the single-MS-antenna case in (4). Consequently, the algorithm proposed for the single-MSantenna case such as DCOMP can be used without any modification in the multiple-MS-antenna
case regardless of different approaches.
B. Extension to wideband systems
In this subsection, we extend our work in the narrowband case into the wideband case. Since
the same algorithm can be used for both the single-MS-antenna case and the multiple-MS-antenna

19

case as shown in Section V-A, we focus on the single antenna case for the sake of exposition.
Compared to the narrowband channel model in (1), we adopt the delay-d MIMO channel
model for the wideband OFDM systems with K subcarriers [42], [43]. Let Ts be the sampling
period, τ` be the delay of the `-th path, NCP be the cyclic prefix length, and p(t) denote a filter
comprising the combined effects of pulse shaping and other analog filtering. The delay-d MIMO
channel matrix is modeled as
ht [d] =

L
X

g`,t p(dTs − τ` )a(φ` ) for d = 0, ..., NCP−1 ,

(38)

`=1

and the channel frequency response matrix at each subcarrier k can be expressed as
ht,k =

L
X

g`,t c`,k a(φ` ),

(39)

`=1

where c`,k =

PNCP −1
d=0

p(dTs − τ` )e−

j2πkd
K

. This can be cast in the compressive sensing framework

as ht,k = A (gt } ck ), where gt and ck are sparse vectors with L nonzero elements associated
with g`,t and c`,k . Note that gt and ck share the same support for all t and k.
An analog combiner WRF,t at frame t must be common for all subcarrier k. Since Proposition 2
shows that changing baseband combiners does not impact the performance as long as the analog
combiner is fixed, we use a common sensing matrix Wt = WBB,k WRF,t for all subcarriers at
frame t. Then, the signal part of the received signal in the baseband at frame t and subcarrier
k is given by rt,k = Wt A (gt } ck ).
The problem of finding gt } ck is similar to the narrowband case in Section III because
gt } ck share the same support for all t and k. Therefore, the DCOMP algorithm in Section III
can be directly applied to this problem. This fails, however, to exploit a useful property of sparse
wideband channels. At a fixed subcarrier k = k0 , the time domain average of rt,k0 r∗t,k0 becomes
T
X

rt,k0 r∗t,k0

=

t=1

T
X

Wt A (gt } ck0 ) (gt } ck0 )∗ A∗ Wt∗ .

(40)

t=1

Since (40) has a similar form to (17) of the narrowband case, DCOMP can be used for this time
domain operation. The frequency domain average of rt0 ,k r∗t0 ,k at a fixed frame t = t0 can be
expressed differently from (40) as
K
X
k=1

rt0 ,k r∗t0 ,k = Wt0 A

K
X
k=1

!
(gt0 } ck ) (gt0 } ck )∗ A∗ Wt∗0 ,

(41)

20

1

1

0.9

0.9

0.8

0.8

0.7

0.7
0.6

T = 1, 4, 8, 64, ∞

CDF

CDF

0.6
0.5
0.4

0.5
0.4

0.3

Lo = 0, 3, 7
(SOMP)

0.3
T = 1 (OMP)
T =4
T =8
T = 64
Asymptotic (T → ∞)

0.2
0.1

Lo = 0, 3, 7
(OMP)

0.1

0

0
0

0.5

1

1.5

T
T
T
T
T
T

0.2

2

0

0.5

ρ(S) (S, S0 , Φ, G)

1

1.5

2

= 1 (OMP), Lo = 0
= 1 (OMP), Lo = 3
= 1 (OMP), Lo = 7
→ ∞, Lo = 0
→ ∞, Lo = 3
→ ∞, Lo = 7

2.5

3

ρ(S) (S, S0 , Φ, G)

(a) CDF of ρ(S) (S, S0 , Φ, G) when Lo = 0

(b) CDF of ρ(S) (S, S0 , Φ, G) when T = 1 or T = ∞

Fig. 3. CDFs of ρ(S) (S, S0 , Φ, G). The optimal support set S has L = 8 elements, the optimally preselected subset So has Lo
elements, Φ ∈ CM ×N is a fixed sensing matrix of M = 8 and N = 64, and G ∈ CL×T is a random matrix whose elements
have zero-mean and unit-variance.

which is similar to (16) of the narrowband case. Since COMP is developed based on the signal
model in (16), COMP can be used for the frequency domain operation in the wideband case.
The different features between the time and frequency domain motivate us to apply different
approaches to each domain, i.e., COMP to the frequency domain and DCOMP to the time domain.
This combination of two algorithms, which we call WB-DCOMP, is described in Algorithm 6.
It is worthwhile to compare WB-DCOMP to the direct extension of COMP and DCOMP. If
the sensing matrix Φ is varying over both time and frequency, the direct extension of DCOMP
provides the best performance. However, if Φ is varying only in the time domain and fixed over
frequency, WB-DCOMP, which is the combination of DCOMP and COMP, can outperform the
direct extension of each algorithm in the wideband case. This is demonstrated in Section VI.
VI. S IMULATION RESULTS
In this section, we validate our analysis in Section IV in terms of the recovery success
probability in the noiseless case especially when M has a similar value to L. We also present
simulation results to demonstrate the performance of the proposed spatial channel covariance
estimation algorithms for hybrid architecture.
A. Theoretical anaylsis in the noiseless case
We compare OMP and SOMP in terms of the cumulative distribution function (CDF) of the
newly defined metric ρ(S) (S, S0 , Φ, G) in (21). Fig. 3 shows how the CDF of ρ(S) (S, S0 , Φ, G)

1

1

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

CDF

CDF

21

0.5
0.4

0.5
0.4

T = 1 (OMP)
T =4
T =8
T = 64
T = 1024
Upper bound as T → ∞

0.3
0.2
0.1

T = 1 (OMP)
T =4
T =8
T = 64
T = 1024
Upper bound as T → ∞

0.3
0.2
0.1

0

0
0

0.5

1

ρ(DS) (S, S0 , Φ1 , · · · , ΦT , G)

(a) Lo = 0

1.5

2

0

1

2

3

4

5

ρ(DS) (S, S0 , Φ1 , · · · , ΦT , G)

(b) Lo = 7

Fig. 4. CDFs of ρ(DS) (S, S0 , Φ1 , ..., ΦT , G). The optimal support set S has L = 8 elements, the optimally preselected subset
So has Lo elements, Φ1 , ..., ΦT ∈ CM ×N are time-varying sensing matrices of M = 8 and N = 64, and G ∈ CL×T is a
random matrix whose elements have zero-mean and unit-variance. The upper bound of ρ(DS) (S, S0 , Φ1 , ..., ΦT , G) as T → ∞
(DS)
denotes ρT →∞ (N, M, L, Lo ) in Theorem 2.

changes according to T and Lo . In Fig. 3(a), Pr(ρ(S) (S, S0 , Φ, G) < 1) becomes higher as T
(S)

gets larger. In addition, as shown in Proposition 1, the CDF converges to that of ρT →∞ (S, S0 , Φ)
in (21) regardless of G. Note that T = 1 indicates the OMP case. Fig. 3(b) shows the CDF
versus Lo in the OMP case (T = 1) and the asymptotic SOMP case (T = ∞). As Lo increases,
the gap between OMP and SOMP reduces and becomes zero if M = L. In this extreme case
where M = L, both OMP and SOMP do not work properly because Pr(ρ(S) (S, S0 , Φ, G) < 1)
becomes low in the last iteration.
Fig. 4 shows the CDF of ρ(DS) (S, S0 , Φ1 , ..., ΦT , G) in the DSOMP case. Unlike the SOMP
(DS)

case, ρ(DS) (S, S0 , Φ1 , ..., ΦT , G) converges to a deterministic value ρT →∞ (N, M, L, Lo ) in (24)
for any sensing matrices Φ1 , ..., ΦT as T increases. The simulation results in Fig. 4 also indicate
(DS)

that ρT →∞ (N, M, L, Lo ) has an upper bound as shown in (24), which is consistent with the
analytical result in Theorem 2. Although the convergence rate becomes slower as Lo increases,
Pr(ρ(DS) (S, S0 , Φ1 , ..., ΦT , G) < 1) = 1 is guaranteed as T → ∞. Simulation results in Fig. 4
show that a reasonably small value of T  ∞ can ensure a perfect recovery.
Fig. 5 compares DSOMP and DCOMP. As discussed in Section IV, both select the same
support at the first iteration, and the superiority of DCOMP over DSOMP starts from the second
iteration. In Fig. 5, the CDF of ρ(DS) (S, S0 , Φ1 , ..., ΦT , G) and ρ(DC) (S, S0 , Φ1 , ..., ΦT , G) are
compared in the last iteration, i.e., Lo = 7. The figure shows that the CDF of the DCOMP
case is located on the left side of that of the DSOMP case, which means DCOMP has a higher

22

1
0.9
0.8
0.7

CDF

0.6
0.5
T
T
T
T
T
T
T
T

0.4
0.3
0.2
0.1

= 1: DSOMP
= 1: DCOMP
= 4: DSOMP
= 4: DCOMP
= 8: DSOMP
= 8: DCOMP
= 64: DSOMP
= 64: DCOMP

0
0

1

2

3

4

5

ρ(DS or DC) (S, S0 , Φ1 , · · · , ΦT , G)

1

1

0.9

0.9

0.8
0.7
0.6
0.5

1st iter.
2nd iter.
3rd iter.
4th iter.
5th iter.
6th iter.
7th iter.
8th iter.
total

0.4
0.3
0.2
0.1

Success probability per iteration

Success probability per iteration

Fig. 5. Comparison between ρ(DS) (S, S0 , Φ1 , ..., ΦT , G) of DSOMP and ρ(DC) (S, S0 , Φ1 , ..., ΦT , G) of DCOMP when
N = 64, M = 8, L = 8, and Lo = 7.

0.8
0.7
0.6
0.5

1st iter.
2nd iter.
3rd iter.
4th iter.
5th iter.
6th iter.
7th iter.
8th iter.
total

0.4
0.3
0.2
0.1

0

0
0

10

20

30

40

50

60

0

10

T

20

30

40

50

60

T

(a) DSOMP

(b) DCOMP

Fig. 6. Success probability of recovering one of elements in the optimal support set at each iteration for Lo = 1, ..., L when
N = 64, M = 8, and L = 8.

success probability than DSOMP.
Simulation results with respect to the success probability per iteration is shown in Fig. 6. As
expected, the success probability at the first iteration is identical for both DSOMP and DCOMP,
but the success probability of DCOMP becomes higher than that of DSOMP from the second
iteration, leading to the higher success probability in total.
B. Performance evaluation on the covariance estimation for hybrid architecture
In this subsection, we evaluate the performance
of the spatial
channel covariance estimation.


The performance metric is defined as η =

Tr U∗ Rh UR̂
R̂h
h


Tr U∗R Rh URh

[24] where URh and UR̂h are the

h

matrices whose columns are the eigenvectors of the ideal covariance Rh and the estimated
covariance R̂h . Note that 0 ≤ η ≤ 1 and a larger η indicates a more accurate estimation.

1

1

0.9

0.9

0.8

0.8

0.7

0.7

Efficiency metric

Efficiency metric

23

0.6
0.5
0.4
0.3
OMP
SOMP
COMP
DSOMP
DCOMP

0.2
0.1

0.6
0.5
0.4
0.3
OMP
SOMP
COMP
DSOMP
DCOMP

0.2
0.1

0

0
0

10

20

30

40

50

0

T

(a) M = 16

10

20

30

40

50

T

(b) M = 8

Fig. 7. Comparison among OMP, SOMP, COMP, DSOMP, and DCOMP when N = 64, D = 256, L = 8 for the single-MSantenna case in narrowband systems.

Fig. 7 shows simulation results when N = 64, D = 256, L = 8, and SNR = 10 dB for the
single-MS-antenna case in narrowband systems. In the fixed combining matrix case in Fig. 7(a)
where M = 16, we can see that the proposed COMP outperforms OMP and SOMP. Combined
with the time-varying analog combining matrix, DSOMP and DCOMP have more gain over other
techniques with a fixed combining matrix. The gap between DSOMP and DCOMP, however, is
marginal in this case.
Fig. 7(b) shows the results with 8 RF chains instead of 16 RF chains. As discussed in
Section IV, none of the algorithms using a fixed analog combining matrix such as OMP, SOMP,
and COMP work properly. Moreover, SOMP even yields worse performance than OMP. In
contrast, DSOMP and DCOMP that use a time-varying sensing matrix have a remarkable gain
compared to those that use a fixed sensing matrix. In addition, DCOMP, which exploits the
Hermitian property of covariance matrices, has a considerable gain compared to DSOMP. These
results are consistent with the analysis in Section IV.
Regarding the extension of the proposed work to the multiple-MS-antenna case, the four
different approaches explained in Section V-A are compared in Fig. 8. As expected, the use
of time-varying precoding and combining matrices at both BS and MS improves covariance
estimation. Fig. 9 shows the extension to the wideband case. As shown in the figure, the
combination of COMP and DCOMP outperforms the direct extension of COMP or DCOMP. In
addition to the efficiency metric shown in Fig. 8(a) and Fig. 9(a), we evaluate the loss caused
by covariance estimation error in terms of spectral efficiency of hybrid precoding. The rate loss

1

0

0.9

-10

0.8

-20

0.7

-30

Rate loss (%)

Efficiency metric

24

0.6
0.5
0.4

-40
-50
-60

0.3
DCOMP
DCOMP
DCOMP
DCOMP

0.2
0.1

w/
w/
w/
w/

Approach
Approach
Approach
Approach

1
2
3
4

(W:
(W:
(W:
(W:

fixed, f :
fixed, f :
varying,
varying,

-70

fixed)
varying)
f : fixed)
f : varying)

DCOMP
DCOMP
DCOMP
DCOMP

-80

0

w/
w/
w/
w/

Approach
Approach
Approach
Approach

1
2
3
4

(W:
(W:
(W:
(W:

fixed, f :
fixed, f :
varying,
varying,

fixed)
varying)
f : fixed)
f : varying)

-90
0

5

10

15

20

0

5

10

T

15

20

T

(a) Efficiency metric vs. T

(b) Rate loss vs. T

Fig. 8. Comparison among four different approaches when NT = 64, MT = 8, DT = 256, NR = 64, MR = 8, DR = 256,
L = 8, and SNR = 0 dB for the multiple-MS-antenna case in narrowband systems.

1

0

0.9

-5

0.8

Rate loss (%)

Efficiency metric

-10
0.7
0.6
0.5
0.4

-15
-20
-25

0.3
-30
direct extension of COMP
direct extension of DSOMP
direct extension of DCOMP
modified DCOMP (WB-DCOMP)

0.2
0.1

direct extension of COMP
direct extension of DSOMP
direct extension of DCOMP
modified DCOMP (WB-DCOMP)

-35

0

-40
0

5

10

15

20

25

30

0

T

5

10

15

20

25

30

T

(a) Efficiency metric vs. T

(b) Rate loss vs. T

Fig. 9. Comparison among the WB-DCOMP and the direct extensions of other algorithms when N = 64, M = 8, D = 256,
L = 8, K = 128, NCP = 32, and SNR = 0 dB for the single-MS-antenna case in wideband systems. The normalized delays
τ` /Ts ’s are uniformly distributed in (0, NCP ), and p(t) = sinc(t/Ts ) is used assuming that ideal low pass filters are employed.

in Fig. 8(b) and Fig. 9(b) is defined as

SEest −SEideal
SEideal

in percentage where SEest and SEideal are

the spectral efficiency of the estimated and ideal covariance case under the assumption that the
analog precoding matrix is composed of the dominant eigenvectors of the estimated or ideal
spatial channel covariance matrix. The figures show that the trend in the rate loss is consistent
with that in the efficiency metric although the rate does not only depend on the efficiency metric
but also on other factors such as the type of MIMO techniques and the number of streams.
VII. C ONCLUSIONS
In this paper, we proposed spatial channel covariance estimation techniques for the hybrid
MIMO structure. Based on compressive sensing techniques that leverage the spatially sparse

25

property of the channel, we developed covariance estimation algorithms that are featured by two
key ideas. One is to apply a time-varying sensing matrix, and the other is to exploit the Hermitian
property of the covariance matrix. Simulation results showed that the proposed greedy algorithms
outperform prior work, and the benefit of adopting the two ideas becomes more significant as
the number of RF chains becomes smaller. We also analyzed the performance of the proposed
algorithms in terms of recovery success probability. The theoretical analysis indicated that the
success probability approaches one as the number of snapshots increases if a time-varying sensing
matrix is applied. The analysis also proved that using the structured property of the covariance
matrix improves the estimation performance, which is consistent with the simulation results.
A PPENDIX A
P ROOF OF T HEOREM 1
Let Pt,So = Φt,So Φ†t,So . Then, the support selection criterion at the (Lo + 1)-th iteration in the
DSOMP shown in Algorithm 3 becomes
j

(opt)

= arg max
j∈N

T
X

φ∗t,j (I − Pt,So ) Φt,S gt,S

2

t=1

(a)

= arg max

j∈N \So

(b)

= arg max

j∈N \So

T
X

φ∗t,j (I − Pt,So ) Φt,S\So gt,S\So

2

(42)

t=1
T
X
X

2

ψ ∗t,j ψ t,i gt,i ,

t=1 i∈S\So

where (a) comes from Φt,S gt,S = Φt,So gt,So + Φt,S\So gt,S\So and (I − Pt,So ) Φt,So = 0, and (b)
comes from (I−Pt,So )2 = I−Pt,So due to the characteristics of the orthogonal projection matrix
Pt,So . Consequently, one of the elements in the optimal support S is selected, i.e., j (opt) ∈ S \ So
at iteration Lo + 1 if and only if (20) is satisfied.

26

A PPENDIX B
P ROOF OF P ROPOSITION 1
If Φ1 = · · · = ΦT = Φ, ρ(DS) (S, So , Φ1 , ..., ΦT , G) in (20) becomes
 P

maxj∈N \S ψ ∗j ΨS\So T1 Tt=1 gt gt∗ Ψ∗S\So ψ j
 P

ρ(DS) (S, So , Φ, ..., Φ, G) =
maxj∈S\So ψ ∗j ΨS\So T1 Tt=1 gt gt∗ Ψ∗S\So ψ j
(a)

→

maxj∈N \S ψ ∗j ΨS\So Ψ∗S\So ψ j

(43)

maxj∈S\So ψ ∗j ΨS\So Ψ∗S\So ψ j
P
2
maxj∈N \S i∈S\So ψ ∗j ψ i
=
P
2,
maxj∈S\So i∈S\So ψ ∗j ψ i

where (a) comes from the fact that gt,i are independent random variables with zero mean and
unit variance.
A PPENDIX C
P ROOF OF T HEOREM 2
ρ(DS) (S, So , Φ1 , ..., ΦT , G) in (20) can be written as

P P
2
∗
maxj∈N \S T1 Tt=1
ψ
ψ
+
X
t,j
t,j t,i
i∈S\So

P
ρ(DS) (S, So , Φ1 , ..., ΦT , G) =
P
2
∗
ψ
ψ
+
X
maxj∈S\So T1 Tt=1
t,j
t,j t,i
i∈S\So
h
i
2
(L − Lo )E ψ ∗t,j(6=i) ψ t,i
(a)
h
i
h
i,
→
2
2
(L − Lo − 1)E ψ ∗t,j(6=i) ψ t,i + E ψ ∗t,i ψ t,i
where Xt,j =

P

i1 ∈S\So

P

i2 ∈S\So
i2 6=i1

(44)

∗
gt,i1 gt,i
ψ ∗t,j ψ t,i1 ψ ∗t,i2 ψ t,j . In (44), (a) comes from the fact that
2

ψ ∗t,i ψ t,i for all i are identical random variables with non-zero mean, ψ ∗t,j ψ t,i for all j 6= i are
identical random variables with zero-mean, and gt,i are independent random variables with zero
mean and unit variance.
h
i
h
i
2
2
Now let us look at E ψ ∗t,i ψ t,i
and E ψ ∗t,j(6=i) ψ t,i . Let Ωt = Ψ∗t,N \So Ψt,N \So , then the
trace of Ωt becomes
Tr(Ωt ) = Tr Φ∗t,N \So (IM − Pt,So )2 Φt,N \So
(a)



= Tr (IM − Pt,So ) N IM − Φt,So Φ∗t,So

(b)

= N (M − Lo ) ,



(45)

27

where (a) comes from Φt,So Φ∗t,So + Φt,N \So Φ∗t,N \So = N IN due to the tight frame property, and
i
h
2
∗
is given by
(b) comes from (I − Pt,So ) Φt,So = 0. From (45), the lower bound of E ψ t,i ψ t,i
h
E

2
ψ ∗t,i ψ t,i


2
i (a)
2
 ∗
E [Tr(Ωt )]
N 2 (M − Lo )2
=
≥ E ψ t,i ψ t,i
=
,
N − Lo
(N − Lo )2

(46)

where (a) comes from the fact that E[|X|2 ] ≥ (E[|X|])2 for any random variable X.
Now, let us look at the squared Frobenius norm of Ωt that is given by

2 
2
∗
kΩt kF = Tr (IM − Pt,So ) Φt,N \So Φt,N \So

(47)

2

= N (M − Lo ) .
In addition, E [kΩt k2F ] can be represented differently as

 X


E kΩt k2F = E 


2

X

ψ ∗t,i ψ t,i +

i∈N \So

= (N − Lo )E

j,i∈N \So
j6=i

h

2
ψ ∗t,i ψ t,i

i



2
ψ ∗t,j ψ t,i 


(48)

+ (N − Lo )(N − Lo − 1)E

h

2
ψ ∗t,j(6=i) ψ t,i

i

.

h
i
2
From (46)-(48), the upper bound of E ψ ∗t,j(6=i) ψ t,i
can be obtained as
h

E ψ ∗t,j(6=i) ψ t,i

2

i

=

h
i
2
N 2 (M − Lo ) − (N − Lo )E ψ ∗t,i ψ t,i

(N − Lo )(N − Lo − 1)
(49)
N (M − Lo )(N − M )
≤
.
(N − Lo )2 (N − Lo − 1)
i
i
h
h
2
2
By using the lower bound of E ψ ∗t,i ψ t,i
and the upper bound of E ψ ∗t,j(6=i) ψ t,i , it can
2

be shown that the converged value of ρ(DS) (S, So , Φ1 , ..., ΦT , G) in (44) has an upper bound as
ρ(DS) (S, So , Φ1 , ..., ΦT , G) →

L − Lo
L − Lo − 1 +

≤

h
i
2
E |ψ ∗t,i ψ t,i |
h
i
2
E |ψ ∗t,j(6=i) ψ t,i |

L − Lo
L − Lo − 1 +

(M −Lo )(N −Lo −1)
(N −M )

(50)
,

and this upper bound is always less than or equal to one because N ≥ M ≥ L > Lo .

28

A PPENDIX D
P ROOF OF T HEOREM 4
As T → ∞, ρ(DS) (S, So , Φ1 , ..., ΦT , G) in (20) can be rewritten as
ρ

(DS)

(S, So , Φ1 , ..., ΦT , G) =

(DS)
∗
t=1 φt,j Qt,S,So φt,j
P
(DS)
maxj∈S\So Tt=1 φ∗t,j Qt,S,So φt,j

maxj∈N \S

PT

i
h
(DS)
∗
Tr(Φt,N \S Q̃t,S,So Φt,N \S )
(a)
i,
h
→
(DS)
1
∗
Φ
)
Q̃
E
Tr(Φ
t,S\S
o
t,S\So t,S,So
L−Lo
1
E
N −L

(51)

(DS)

where Q̃t,S,So = (IM − Pt,So ) Φt,S\So Φ∗t,S\So (IM − Pt,So ) and (a) comes from the fact that G
and Φt are independent and all elements in G have zero mean and unit variance.
In a similar way, ρ(DC) (S, So , Φ1 , ..., ΦT , G) in (25) converges as
i
h
(DC)
1
∗
E Tr(Φt,N \S Q̃t,S,So Φt,N \S )
N −L
i,
h
ρ(DC) (S, So , Φ1 , ..., ΦT , G) →
(DC)
1
∗
E Tr(Φt,S\So Q̃t,S,So Φt,S\So )
L−Lo

(52)

(DC)

where Q̃t,S,So = Φ∗t,S\So Φ∗t,S\So − Pt,So Φt,S\So Φ∗t,S\So Pt,So .
Since the terms in the expectation in (51) and (52) are independent random variables with
respect to t, we omit the time slot index t for simplicity. The difference between the inner parts
of the denominators in (51) and (52) becomes





(DC)
(DS)
Tr Φ∗S\So Q̃S,So − Q̃S,So ΦS\So = 2Re Tr Φ∗S\So (I − PSo )ΦS\So Φ∗S\So PSo ΦS\So . (53)
Note that both PSo and I − PSo can be represented as the production of two semi-unitary
matrices, and thus Φ∗S\So PSo ΦS\So and Φ∗S\So (I − PSo ) ΦS\So become positive semidefinite
matrices. Since the trace of the product of two semidefinite matrices is larger than or equal
to zero, (53) becomes
Tr



Φ∗S\So



(DC)
Q̃S,So

−

(DS)
Q̃S,So





ΦS\So ≥ 0

(54)

for any Φ, S, and So . The difference between the inner parts of the numerators in (51) and (52)

29

becomes





(DS)
(DC)
Tr Φ∗N \S Q̃S,So − Q̃S,So ΦN \S = 2Re Tr Φ∗N \S (I − PSo ) ΦS\So Φ∗S\So PSo ΦN \S
(a)

= −2Re Tr PSo ΦS\So Φ∗S\So (I − PSo ) ΦS\So Φ∗S\So




(DC)
(DS)
∗
= −Tr ΦS\So Q̃S,So − Q̃S,So ΦS\So


(55)

≤ 0,
where (a) can be proved by using ΦN \S Φ∗N \S = N IM −ΦSo Φ∗So −ΦS\So Φ∗S\So and PSo (N IM −
ΦSo Φ∗So )(IM − PSo ) = 0.
From the two inequalities in (54) and (55), the converged value in (51) is always larger than
or equal to that in (52), and this completes the proof.
R EFERENCES
[1] W. Roh, J. Seol, J. Park, B. Lee, J. Lee, Y. Kim, J. Cho, K. Cheun, and F. Aryanfar, “Millimeter-wave beamforming as
an enabling technology for 5G cellular communications: theoretical feasibility and prototype results,” IEEE Comm. Mag.,
vol. 52, no. 2, pp. 106–113, Feb. 2014.
[2] R. Heath, N. González-Prelcic, S. Rangan, W. Roh, and A. Sayeed, “An overview of signal processing techniques for
millimeter wave MIMO systems,” IEEE Jour. Select. Topics in Sig. Proc., vol. 10, no. 3, pp. 436–453, Apr. 2016.
[3] T. Marzetta, “Noncooperative cellular wireless with unlimited numbers of base station antennas,” IEEE Trans. Wireless
Comm., vol. 9, no. 11, pp. 3590–3600, Nov. 2010.
[4] J. Hoydis, S. ten Brink, and M. Debbah, “Massive MIMO in the UL/DL of cellular networks: How many antennas do we
need?” IEEE Jour. Select. Areas in Comm., vol. 31, no. 2, pp. 160–171, Feb. 2013.
[5] A. F. Molisch, V. V. Ratnam, S. Han, Z. Li, S. L. H. Nguyen, L. Li, and K. Haneda, “Hybrid beamforming for massive
MIMO: A survey,” IEEE Comm. Mag., vol. 55, no. 9, pp. 134–141, 2017.
[6] O. El Ayach, S. Rajagopal, S. Abu-Surra, Z. Pi, and R. Heath, “Spatially sparse precoding in millimeter wave MIMO
systems,” IEEE Trans. Wireless Comm., vol. 13, no. 3, pp. 1499–1513, Mar. 2014.
[7] A. Alkhateeb, O. El Ayach, G. Leus, and R. Heath, “Hybrid precoding for millimeter wave cellular systems with partial
channel knowledge,” in Proc. Info. Th. App. (ITA), Feb. 2013, pp. 1–5.
[8] X. Yu, J. Shen, J. Zhang, and K. Letaief, “Alternating minimization algorithms for hybrid precoding in millimeter wave
MIMO systems,” IEEE Jour. Select. Topics in Sig. Proc., vol. 10, no. 3, pp. 485–500, Apr. 2016.
[9] S. Park, A. Alkhateeb, and R. Heath, “Dynamic subarrays for hybrid precoding in wideband mmWave MIMO systems,”
IEEE Trans. Wireless Comm., vol. 16, no. 5, pp. 2907–2920, May 2017.
[10] W. Ni and X. Dong, “Hybrid block diagonalization for massive multiuser MIMO systems,” IEEE Trans. Comm., vol. 64,
no. 1, pp. 201–211, Jan. 2016.
[11] F. Sohrabi and W. Yu, “Hybrid digital and analog beamforming design for large-scale antenna arrays,” IEEE Jour. Select.
Topics in Sig. Proc., vol. 10, no. 3, pp. 501–513, Apr. 2016.
[12] A. Alkhateeb, G. Leus, and R. Heath, “Limited feedback hybrid precoding for multi-user millimeter wave systems,” IEEE
Trans. Wireless Comm., vol. 14, no. 11, pp. 6481–6494, Nov. 2015.

30

[13] A. Adhikary, J. Nam, J. Ahn, and G. Caire, “Joint spatial division and multiplexing: The large-scale array regime,” IEEE
Trans. Info. Th., vol. 59, no. 10, pp. 6441–6463, Oct. 2013.
[14] S. Park, J. Park, A. Yazdan, and R. W. Heath, “Exploiting spatial channel covariance for hybrid precoding in massive
MIMO systems,” IEEE Trans. Sig. Proc., vol. 65, no. 14, pp. 3818–3832, Jul. 2017.
[15] A. Alkhateeb, O. El Ayach, G. Leus, and R. Heath, “Channel estimation and hybrid precoding for millimeter wave cellular
systems,” IEEE Jour. Select. Topics in Sig. Proc., vol. 8, no. 5, pp. 831–846, Oct. 2014.
[16] J. Lee, G. T. Gil, and Y. H. Lee, “Channel estimation via orthogonal matching pursuit for hybrid MIMO systems in
millimeter wave communications,” IEEE Trans. Comm., vol. 64, no. 6, pp. 2370–2386, Jun. 2016.
[17] R. T. Suryaprakash, M. Pajovic, K. J. Kim, and P. Orlik, “Millimeter wave communications channel estimation via Bayesian
group sparse recovery,” in Proc. IEEE Int. Conf. Acoust., Speech and Sig. Proc., Mar. 2016, pp. 3406–3410.
[18] Z. Li, S. Han, and A. F. Molisch, “Optimizing channel-statistics-based analog beamforming for millimeter-wave multi-user
massive MIMO downlink,” IEEE Trans. Wireless Comm., vol. 16, no. 7, pp. 4288–4303, Jul. 2017.
[19] R. Méndez-Rial, N. G. Prelcic, and R. W. Heath, “Adaptive hybrid precoding and combining in mmWave multiuser MIMO
systems based on compressed covariance estimation,” in Proc. IEEE Int. Workshop on Comp. Adv. in Multi-Sensor Adap.
Proc. (CAMSAP), Dec. 2015, pp. 213–216.
[20] D. Love, R. Heath, V. Lau, D. Gesbert, B. Rao, and M. Andrews, “An overview of limited feedback in wireless
communication systems,” IEEE Jour. Select. Areas in Comm., vol. 26, no. 8, pp. 1341–1365, Oct. 2008.
[21] E. Bjornson, D. Hammarwall, and B. Ottersten, “Exploiting quantized channel norm feedback through conditional statistics
in arbitrarily correlated MIMO systems,” IEEE Trans. Sig. Proc., vol. 57, no. 10, pp. 4027–4041, Oct. 2009.
[22] S. Shakeri, D. D. Ariananda, and G. Leus, “Direction of arrival estimation using sparse ruler array design,” in Proc. IEEE
Workshop on Sign. Proc. Adv. in Wireless Comm., June 2012, pp. 525–529.
[23] D. Romero and G. Leus, “Compressive covariance sampling,” in Proc. Info. Th. App. (ITA), Feb 2013, pp. 1–8.
[24] S. Haghighatshoar and G. Caire, “Massive MIMO channel subspace estimation from low-dimensional projections,” IEEE
Trans. Sig. Proc., vol. 65, no. 2, pp. 303–318, Jan 2017.
[25] H. Ghauch, T. Kim, M. Bengtsson, and M. Skoglund, “Subspace estimation and decomposition for large millimeter-wave
MIMO systems,” IEEE Jour. Select. Topics in Sig. Proc., vol. 10, no. 3, pp. 528–542, Apr. 2016.
[26] Y. Chi, Y. C. Eldar, and R. Calderbank, “PETRELS: parallel subspace estimation and tracking by recursive least squares
from partial observations,” IEEE Trans. Sig. Proc., vol. 61, no. 23, pp. 5947–5959, Dec 2013.
[27] S. F. Chuang, W. R. Wu, and Y. T. Liu, “High-resolution AoA estimation for hybrid antenna arrays,” IEEE Trans. on Ant.
Prop., vol. 63, no. 7, pp. 2955–2968, July 2015.
[28] Y. Peng, Y. Li, and P. Wang, “An enhanced channel estimation method for millimeter wave systems with massive antenna
arrays,” IEEE Comm. Lett., vol. 19, no. 9, pp. 1592–1595, Sep. 2015.
[29] J. Lee, G. T. Gil, and Y. H. Lee, “Channel estimation via orthogonal matching pursuit for hybrid MIMO systems in
millimeter wave communications,” IEEE Trans. Comm., vol. 64, no. 6, pp. 2370–2386, Jun. 2016.
[30] S. F. Cotter, B. D. Rao, K. Engan, and K. Kreutz-Delgado, “Sparse solutions to linear inverse problems with multiple
measurement vectors,” IEEE Trans. Sig. Proc., vol. 53, no. 7, pp. 2477–2488, July 2005.
[31] J. Chen and X. Huo, “Theoretical results on sparse representations of multiple-measurement vectors,” IEEE Trans. Sig.
Proc., vol. 54, no. 12, pp. 4634–4643, Dec. 2006.
[32] M. F. Duarte and Y. C. Eldar, “Structured compressed sensing: From theory to applications,” IEEE Trans. Sig. Proc.,
vol. 59, no. 9, pp. 4053–4085, Sep. 2011.
[33] J. F. Determe, J. Louveaux, L. Jacques, and F. Horlin, “On the exact recovery condition of simultaneous orthogonal
matching pursuit,” IEEE Sig. Proc. Lett., vol. 23, no. 1, pp. 164–168, Jan 2016.

31

[34] Y. C. Eldar and G. Kutyniok, Compressed Sensing: Theory and Applications.

Cambridge, 2012.

[35] S. Park and R. W. Heath, “Spatial channel covariance estimation for mmWave hybrid MIMO architecture,” in Proc. of
Asilomar Conf. on Sign., Syst. and Computers, Nov. 2016, pp. 1424–1428.
[36] R. Méndez-Rial, C. Rusu, N. G. Prelcic, A. Alkhateeb, and R. W. Heath, “Hybrid MIMO architectures for millimeter wave
communications: Phase shifters or switches?” IEEE Access, vol. 4, pp. 247–267, Jan. 2016.
[37] J. A. Tropp, A. C. Gilbert, and M. J. Strauss, “Simultaneous sparse approximation via greedy pursuit,” in Proc. IEEE Int.
Conf. Acoust., Speech and Sig. Proc., vol. 5, Mar. 2005, pp. v/721–v/724 Vol. 5.
[38] J. Determe, J. Louveaux, L. Jacques, and F. Horlin, “On the noise robustness of simultaneous orthogonal matching pursuit,”
IEEE Trans. Sig. Proc., vol. 65, no. 4, pp. 864–875, Feb. 2017.
[39] M. A. Davenport and M. B. Wakin, “Analysis of orthogonal matching pursuit using the restricted isometry property,” IEEE
Trans. Info. Th., vol. 56, no. 9, pp. 4395–4401, Sept 2010.
[40] D. D. Ariananda and G. Leus, “Direction of arrival estimation for more correlated sources than active sensors,” Sig. Proc.,
vol. 93, no. 12, pp. 3435–3448, 2013.
[41] J. A. Tropp, “Greed is good: algorithmic results for sparse approximation,” IEEE Trans. Info. Th., vol. 50, no. 10, pp.
2231–2242, Oct 2004.
[42] P. Schniter and A. Sayeed, “Channel estimation and precoder design for millimeter-wave communications: The sparse
way,” in Proc. of Asilomar Conf. on Sign., Syst. and Computers, Nov. 2014, pp. 273–277.
[43] A. Alkhateeb and R. Heath, “Frequency selective hybrid precoding for limited feedback millimeter wave systems,” IEEE
Trans. Comm., vol. 64, no. 5, pp. 1801–1818, May 2016.

