1

Low-Complexity Massive MIMO Subspace
Estimation and Tracking from Low-Dimensional
Projections

arXiv:1608.02477v3 [] 24 Jul 2017

Saeid Haghighatshoar, Member, IEEE, Giuseppe Caire, Fellow, IEEE

Abstract—Massive MIMO is a variant of multiuser MIMO,
in which the number of antennas M at the base-station is
very large, and generally much larger than the number of
spatially multiplexed data streams to/from the users. It has
been observed that in many realistic propagation scenarios
as well as in spatially correlated channel models used in
standardizations, although the user channel vectors have a
very high-dim M , they lie on low-dim subspaces due to
their limited angular spread (spatial correlation). This low-dim
subspace structure remains stable across many coherence blocks
and can be exploited in several ways to improve the system
performance. A main challenge, however, is to estimate this
signal subspace from samples of users’ channel vectors as fast
and efficiently as possible. In a recent work, we addressed this
problem and proposed a very effective novel algorithm referred
to as Approximate Maximum-Likelihood (AML), which was
formulated as a semi-definite program (SDP). In this paper,
we address two problems left open in our previous work,
namely, computational complexity and tracking. The algorithm
proposed in this paper is reminiscent of Multiple Measurement
Vectors (MMV) problem in Compressed Sensing and is proved
to be equivalent to the AML Algorithm for sufficiently dense
angular grids. It has also a very low computational complexity
and is able to track sharp transitions in the channel statistics
very quickly. Although mainly motivated by massive MIMO
applications, our proposed algorithm is of independent interest
in other related subspace estimation applications. We assess
the estimation/tracking performance of our proposed algorithm
empirically via numerical simulations, especially in practically
relevant situations where a direct implementation of the
SDP would be infeasible in real-time. We also compare the
performance of our algorithm with other related subspace
estimation/tracking algorithms in the literature.

I. I NTRODUCTION
ONSIDER a multiuser massive MIMO system formed
by a base-station (BS) with M antennas serving K
single-antenna mobile users in a cellular system. Following
the current massive MIMO approach [1–4], we focus on uplink
(UL) and downlink (DL) in Time Division Duplexing (TDD),
where the base-station (BS) transmit/receive hardware is
designed or calibrated in order to preserve UL-DL reciprocity
[5, 6] such that the BS estimates the channel vectors of the
users from UL orthogonal training pilots sent by the users and
uses them to transmit data to the users in the DL via coherent
beamforming. Since there is no multiuser interference in the

C

The authors are with the Communications and Information Theory Group,
Technische Universität Berlin ({saeid.haghighatshoar, caire}@tu-berlin.de).
A shorter version of this paper was presented in the IEEE International
Conference on Communications (ICC), Paris, France, May 2017.

UL training phase (after neglecting the pilot contamination),
in this paper we focus on the basic channel estimation problem
for a single user.
In massive MIMO systems, the number of antennas M ,
thus, the dimension of the received signal at the BS is very
large. However, in many relevant scenarios, channel vectors
of each user are spatially correlated since the propagation
between the user and the BS occurs through a small set
of Angle of Arrivals (AoAs). Denoting by h(t) ∈ CM the
channel vector of a generic user at a time slot t, this implies
that the signal covariance matrix of h(t) given by S :=
E[h(t)h(t)H ] is typically low-rank. This low-rank structure
can be exploited to improve the system multiplexing gain
and decrease the training overhead. A particularly effective
scheme is the Joint Spatial Division and Multiplexing (JSDM)
approach proposed and analyzed in [7–9], where the users
are partitioned into G > 1 groups such that users in each
group have similar channel subspaces [7–9]. These groups are
separated by a zero-forcing beamforming that uses only the
group subspace information and reduces the dimensionality for
each group g to some mg  M . Then, additional multiuser
multiplexing gain in each group g is obtained by applying
the conventional linear precoding to the lower-dim projected
channel. This has the additional non-trivial advantage that only
m  M RF chains (A/D converters and modulators) are
needed, thus, reducing the A/D conversion rate significantly.
This and many other related examples evidently illustrate
that estimation of the signal subspace of the users plays a
crucial rule in massive MIMO systems. However, obtaining the
signal subspace information in massive MIMO is a high-dim
estimation problem and becomes quite challenging as the
number of antennas M increases, especially that, due to the
limited number of available RF chains at the receiver front end,
the subspace estimation needs to be done with only low-dim
projections of the received signal.
A. Contribution
In our recent work [10, 11], we studied this problem and
developed a new family of efficient algorithms for subspace
estimation in massive MIMO. We also demonstrated via
numerical simulations that our proposed algorithms provide
near-ideal performance for a massive MIMO JSDM system.
However, the low-complexity implementation of our proposed
algorithms was left open in [10, 11]. In this paper, we bridge
the complexity gap by providing efficient and low-complexity
implementation of the algorithms in our previous work

2

[10, 11], with a special focus on the AML (approximate
maximum likelihood) Algorithm. Our approach is based
on approximating the typically high-complexity semi-definite
program (SDP) proposed for the original form of AML
Algorithm in [10, 11] with another convex optimization
problem that can be efficiently solved. We consider a
generalization of the originally proposed AML Algorithm
where the projection (sampling) operator may be time-variant,
i.e., changing in different training slots. This results in
further improvement in the subspace estimation. We extend
our proposed low-complexity algorithm to more practical
array configurations such as 2D rectangular lattice arrays,
and provide guidelines for efficient numerical implementation
for general array configurations. We also illustrate that our
algorithm can be run in a tracking mode, where the subspace
estimate is updated upon arrival of a new training sample.
B. Related Work
Subspace estimation and tracking arises in a variety
of problems in signal processing. In general, whenever a
high-dim signal is generated by a linear process that is
governed by a small number of parameters, it can be
represented by a low-dim structure embedded in a higher-dim
space. This occurs in a wide range of applications such
as Direction-of-Arrival (DoA) estimation [12–14], source
localization [15], anomaly detection [16], adaptive filtering
[17], and wireless communication [18]. The dominant signal
subspace is obtained by computing the covariance matrix of
the data (i.e., its second order statistics), e.g., the channel
covariance matrix S = E[h(t)h(t)H ] in massive MIMO
applications, and calculating its Singular Value Decomposition
(SVD), and is classically known as Principal Component
Analysis (PCA) in statistics [19] and Karhunen-Loève
Transform (KLT) in stochastic processes [20]. In practice,
however, the data covariance matrix is not a priori known
and should be estimated from the observed data samples.
In particular, due to complexity reasons, this needs to be
done by taking as few data samples and by consuming
as less storage as possible. This has motivated a vast
line of research on developing efficient and low-complexity
subspace estimation/tracking algorithms [21–25]. Recently by
the advent of Compressed Sensing [26, 27], there has been a
new serge of interest in exploiting low-dim signal structures
such as sparsity and low-rankness, which has revitalized the
popularity of subspace techniques in a variety of problems
including Matrix Completion [28, 29], Super-Resolution [30],
compressive spectral estimation [31, 32], and sparse signal
reconstruction [33–37]
We compare the performance of our proposed algorithm
with the Singular Value Thresholding (SVT) Algorithm [38]
in the batch mode and with PETRELS Algorithm [25] in
the online tracking mode. Both algorithms provide the state
of the art performance in subspace estimation (SVT) and
tracking (PETRELS). Our numerical simulations in Section
VII illustrate that our algorithm performs better that SVT in
the batch mode and is able to track sharp transitions in signal
statistics much faster than PETRELS in the online mode.

Sparse Scattering Channel
(M − 1)d •
4d •

User

3d • θi

..
.

2d •
d•
0•

Fig. 1: Array configuration in a multi-antenna receiver in the
presence of a sparse scattering channel with only few scatterers with
discrete angle of arrivals.

C. Notation
We show vectors by boldface small letters (e.g., x), matrices
by boldface capital letters (e.g., X), scalar constant by
non-boldface letters (e.g., x or X), and sets by calligraphic
letters (e.g., X ). The i-th element of a vector x and the
(i, j)-th element of a matrix X is denoted by [x]i and [X]i,j .
We represent the i-th row and j-th column of a matrix X
with a row vector Xi,. and a column vector X.,j . We denote
the Hermitian and the transpose of a matrix (or a vector)
X by XH and XT . We use tr(.) for the trace operator.
We denote the complex/real inner product of two matrices
(or two vectors) X and Y by hX, Yi = tr(XH Y), and
hX, YiR = Re[hX, Yi]. We use kxk for the l2 -norm of a
vector x, and kXk = hX, Xi = hX, XiR for the Frobenius
norm of a matrix X. We denote a k × k diagonal matrix
with diagonal elements s1 , . . . , sk with diag(s1 , . . . , sk ). We
indicate the output of any optimization problem such as
arg minX∈X f (X) with X∗ . The identity matrix of order p is
denoted by Ip . For an integer k, we use the shorthand notation
[k] for {1, . . . , k}. We denote the big-O notation by O(.).
II. BASIC S ETUP
A. Array and Signal Model
Consider a BS with a large Uniform Linear Array (ULA)
with M  1 antennas. The geometry of the array is shown
in Fig. 1, with antenna elements having a uniform spacing
d. We assume that the BS serves the users in the angular
range [−θmax , θmax ] for some θmax ∈ (0, π/2), and set
λ
d = 2 sin(θ
, where λ = fc0c denotes the wave-length, where
max )
fc is the carrier frequency and c0 is speed of the light. We
consider a simple propagation model in which the transmission
between a user and the BS occurs through p scatterers (see
Fig. 1). The results can be simply extended to a general
scattering model with a general mixed-type (continuous and
discrete) power distribution over the AoA domain as in [11].
One snapshot of the received signal in a window of training
pilots of size T is given by
y(t) = h(t)z(t) + n(t) :=

p
X
l=1

a(θl )wl (t) z(t) + n(t), (1)

3

Pp
where h(t) := l=1 a(θl )wl (t) ∈ CM denotes the channel
vector of the user1 , z(t) ∈ C is the transmitted pilot symbol
of the user, which typically belongs to a signal constellation
such as QAM, wl (t) ∼ CN (0, σl2 ) is the channel gain of
the l-th multipath component, n(t) ∼ CN (0, σ 2 IM ) is the
additive white Gaussian noise of the antenna elements, and
where a(θ) ∈ CM is the array response at AoA θ, whose k-th
component is given by
[a(θ)]k = ejk

2πd sin(θ)
λ

sin(θ)

= ejkπ sin(θmax ) .

(2)

According to the well-known and widely-accepted Wide
Sense Stationary Uncorrelated Scattering (WSSUS) model,
the channel gains of different paths, i.e., {wl (t)}pl=1 , at every
time t ∈ [T ], are uncorrelated [18]. Without loss of generality,
we assume that the pilot symbol z(t) transmitted along the
channel vector h(t) is normalized to z(t) = 1 in all training
snapshots, thus, letting A = [a(θ1 ), . . . , a(θp )], we have
y(t) = h(t) + n(t) = Aw(t) + n(t),

t ∈ [T ],

(3)

where w(t) = (w1 (t), . . . , wp (t))T for different t ∈ [T ]
are statistically independent. Also, we assume that the AoAs
{θl }pl=1 remain invariant over the whole training period of
length T slots. From (3), the covariance of y(t) is given by
Cy = AΣAH + σ 2 IM =

p
X

σl2 a(θl )a(θl )H + σ 2 IM , (4)

l=1

where Σ = diag(σ12 , . . . , σp2 ) is the covariance
matrix of
Pp
2
H
w(t), and where S = E[h(t)h(t)H ] =
l=1 σl a(θl )a(θl )
is the covariance matrix of the channel vectors. It is not
difficult to check that S is a Hermitian positive semi-definite
(PSD) Toeplitz matrix of rank p, where typically p  M . In
practice, the channel vectors are formed by the superposition
of AoA contributions weighted according to a mixed-type
measure γ(dθ) as in [11], containing both discrete masses
in correspondence of specular reflectors and a continuous
measure corresponding
to scattering clusters. In this case,
R
S is given by γ(dθ)a(θ)H a(θ) and is generally full rank
(algebraic rank). However, since γ has a limited angular
support in practice, S can be well approximated by a low-rank
matrix (effective rank), thus, the low-rank assumption is still
valid. The AML Algorithm in [10, 11] and also our proposed
low-complexity scheme in this paper apply to this general case.
B. Sampling Operator
As explained in the introduction, in massive MIMO systems,
it is crucial to be able to recover the signal subspace
1 In this paper, h(t) denotes the channel vector of a generic user at a specific
subcarrier over an OFDM symbol transmitted at a time slot t ∈ [T ] (a specific
time-frequency tile). Although the instantaneous channel vectors, which are
used to transmit/receive data to/from the users in the DL/UL via coherent
beamforming, might be highly correlated, we always assume that the channel
vectors used for subspace estimation are sufficiently separated in time such
that h(t), t ∈ [T ] are independent, where they are also identically distributed
(i.i.d.) due to channel stationarity in time. Even in a fixed time slot t ∈
[T ], one can obtain i.i.d. realizations of the channel vector by sampling at
sufficiently separated subcarriers in the frequency domain since the channel is
also stationary in frequency with a frequency-invariant second order statistics.
Thus, having multiple subcarriers in the frequency domain has the same effect
as taking more observations in the time domain.

of the users from low-dim projections of their received
channel vectors. In general, low-dimensional projections can
be obtained via a m × M matrix B for some m  M , which
can be implemented as part of the analog receiver front-end. A
particularly simple and attractive choice is “antenna selection”,
where B is a binary 0-1 selection matrix with a single
element equal to 1 in each row. In this paper, we always
consider such an antenna selection scheme as the projection
operator, where in each training slot, the BS samples the
output signal of only m  M random antenna elements
via m available RF chains. Also, we consider a general case
in which the antenna selection can be time-variant. Letting
It = {i1 (t), . . . , im (t)} ⊆ [M ] be the indices of m randomly
selected antenna elements at time t ∈ [T ], we denote the
m × M projection matrix by B(t), where the single 1 in
each row is given by [B(t)]k,ik (t) = 1, for k ∈ [m]. Note
that for such a matrix B(t), we have that B(t)B(t)H = Im .
We define the noisy projection (sketch) at time t ∈ [T ] by
e (t), where y(t) is given by
x(t) := B(t)y(t) = B(t)h(t) + n
e (t) := B(t)n(t) ∼ CN (0, σ 2 Im ) denotes the
(3) and where n
projected noise vector.
C. Performance Metric

Our goal is to find an estimate of the dominant signal
subspace of the covariance matrix S of the channel vectors
b be such an estimate and let S = UΛUH
h(t), t ∈ [T ]. Let S
H
b
b
b
b
b We always
and S = UΛU denote the SVD of S and S.
use the convention that the singular values (e.g., those in Λ
b are sorted in a non-increasing order. We define the
and Λ)
normalized power distribution for S by a vector p ∈ RM
+,
where [p]i = PMλi λ with λk denoting the k-th singular
j=1 j
b = [b
b M ] where u
b i denotes the i-th
value of S. Let U
u1 , . . . , u
b We denote the power captured by columns of U
b
column of U.
H
H
2
M
biu
b i i = E[|b
ui h(t)| ]
by a vector q ∈ R+ , where [q]i := hS, u
gives the amount of power of S captured
by
the
1-dim
(rank-1)
PM
bU
b Hi =
biu
bH
projection operator u
.
Note
that
[q]
=
hS,
U
i
i
i=1
2
tr(S) = E[kh(t)k ], which gives the whole power contained
in S. We normalize q and define the estimated normalized
b ∈ RM
power distribution p
p]i = PM[q]i[q] . Let
+ , where [b
j
j=1
Pk
Pk
ηp (k) :=
[p]
and
η
(k)
=
[b
p
]
,
for
k
∈
[M ],
i
i
b
p
i=1
i=1
denote the whole normalized signal power contained in the
b . Note that since U is the
first k component of p and p
SVD basis for S, we always have ηp (k) ≥ ηpb (k), for
b is always
every k ∈ [M ], which implies that the vector p
majorized by p [39]. Also, due to the normalization, we have
ηp (M ) = ηpb (M ) = 1.
In subspace estimation applications in massive MIMO, e.g.,
in JSDM, the goal is to design for each user a low-dim
beamformer that captures a significant amount of the power
of its channel vectors2 . An appropriate distortion measure
η (k)−η (k)
b ) := maxk∈[M ] p ηp (k)pb ,
for such applications is ν(p, p
which captures the maximum ratio of power loss incurred by
beamforming to the dominant k-dim subspace of the estimate
b rather than the optimal k-dim subspace of S, for any
S
b ) := 1 − ν(p, p
b)
arbitrary k ∈ [M ]. We will use Γ(p, p
2 Namely, an M × q matrix V, for some q  M , satisfying VH V = I ,
q




and E kVH h(t)k2 ≥ (1 − )E kh(t)k2 , for some small  ∈ (0, 1).

4

as the metric for assessing the performance of the subspace
b ) ∈ [0, 1], where Γ(p, p
b) = 1
estimation. Note that Γ(p, p
b for some µ > 0. In particular, if
if and only if S = µS
b ) ≥ 1 − , for some fixed  ∈ (0, 1), then the M × k
Γ(p, p
b k ] obtained from the estimate
beamforming matrix [b
u1 , . . . , u
b for an arbitrary k ∈ [M ] is at least (1 − )-optimal with
S
respect to the best k-dim subspace of true signal covariance
matrix S.
III. P ROBLEM S TATEMENT
In this section, we briefly explain the AML Algorithm
proposed for subspace estimation in [10, 11]. For simplicity
of explanation, we assume that the 0-1 sampling operator is a
fixed m × M matrix B for all t ∈ [T ]. We will later consider
the generalized time-variant sampling operator B(t).
Let y(t) = h(t) + n(t) be the noisy user channel vector
received at the array at time t ∈ [T ], and let x(t) = By(t) be
its m-dim projection via B. We assume that the noise variance
σ 2 is known and multiply all the signals by σ1 to normalize the
noise power to 1. For simplicity, we still use
the same notation
H
bx = 1 P
for the normalized signals. Let C
t∈[T ] x(t)x(t) be
T
the sample covariance matrix of the sketches x(t), t ∈ [T ], let
b x = VDVH be its SVD, and define ∆ := C
b 1/2
= VD1/2 .
C
x
The AML Algorithm in [10, 11] is cast as the following SDP:
(S∗ ,K∗ ) =

arg min

tr(BMBH ) + tr(K)

M∈T+ ,K∈Cm×m

subject to



Im + BMBH
∆H

∆
K



 0,

(5)

where T+ denotes the space of all M × M Hermitian PSD
Toeplitz matrices. The optimal solution S∗ of (5) gives an
estimate of the covariance matrix S of the channel vectors.
In [11], we illustrated via numerical simulations that
AML Algorithm has an excellent performance for estimating
user signal subspaces, especially in a JSDM system setup.
Unfortunately, the SDP (5) proposed for AML Algorithm in
[10, 11] is quite time-consuming, especially for a large array
size M . In this paper, instead of directly solving the SDP (5),
we approximate it by another convex optimization problem for
which we provide an efficient and low-complexity algorithm.
Our algorithm works for the more general setup in which the
m × M sampling matrices B(t) may be time-variant, and can
be applied to more practical array configurations such as 2D
rectangular arrays.
IV. M ATHEMATICAL F ORMULATION
A. Equivalent Convex Optimization
Let G be a discrete grid of size G consisting
of the angles

θi := sin−1 (−1 + 2(i−1)
)
sin(θ
)
,
for
i
∈
[G],
over the
max
G
angular range [−θmax , θmax ]. Let G = [a(θ1 ), . . . , a(θG )] be
an M ×G matrix consisting of array responses at AoAs θi ∈ G,
i ∈ [G]. We assume that G is dense enough such that every
signal covariance matrix S can be well approximated by
S ≈ G diag(s1 , . . . , sG )GH =

G
X
i=1

si a(θi )a(θi )H ,

(6)

with appropriate si ≥ 0, i ∈ [G]3 . For a ULA of size M ,
taking G ≈ 2M is typically sufficient for this approximation
to hold. Also, note that any M × M matrix of the form
(6) is a valid channel covariance matrix since it corresponds
to the covariance matrix of the channel vector h(t) =
P
G
i=1 wi (t)a(θi ) consisting of G scatterers with channel gains
wi (t) ∼ CN (0, si ) located at AoA θi .
Now, consider the following convex optimization problem
to be solved for the G × T matrix W:
√
1
(7)
W∗ = arg min kǦW − Xk2 + T kWk2,1 ,
2
W
where kWk2,1 denotes the l2,1 -norm of W defined by
PG
kWk2,1 :=
i=1 kWi,. k, where X = [x(1), . . . , x(T )] is
the m × T matrix of noisy sketches, and where Ǧ = √1m BG
is an m × G matrix with columns of unit l2 -norm. We prove
the following result.
Proposition 1: Assume that the grid G is dense enough
such that every covariance matrix S can be well approximated
according to (6). Then, the SDP (5) and the convex
optimization (7) are equivalent, in the sense that if W∗ is
the minimizer of (7), then the optimal solution of (5) can
be approximated by S∗ = G diag(s∗1 , . . . , s∗G )GH , where
kW∗ k
s∗i = m√i,.T .

Proof: Proof in Appendix A.
Some remarks are in order here.
Remark 1: Optimization problems of the type (7) with
an l2,1 -norm regularization of the form τ kWk2,1 , for some
regularization factor τ > 0, are quite well-known and are
widely applied to solve Multiple Measurement Vector (MMV)
problem in Compressed Sensing [36, 37], where multiple
measurement vectors correspond to different realizations of
a sparse vector all having the same sparsity pattern (location
of nonzero coefficients). It is also well-known that l2,1 -norm
regularization promotes the block sparsity or row sparsity of
the optimal solution W∗ of the channel coefficients in (7). As
∗
each row Wi,.
of W∗ corresponds to the channel gain of a
scatterer located at θi ∈ G across T training slots, considering
the sparse scattering channel in the angular domain, this seems
to be quite a reasonable regularization. However, the main
novelty in (7) consists √
of the remarkable fact that for the
particular choice λ = T of the regularization coefficient
and within the assumptions of Proposition 1, this particular
instance of MMV is, asymptotically for sufficiently dense
angular grids, equivalent to the AML Algorithm, which is
derived in a completely different way without any assumption
on grid quantization in the angular domain.
♦
Remark 2: Proposition 1 implies that by increasing the
number of grid points G the optimal solution S∗ of SDP
(5) can be better approximated as G diag(s1 , . . . , sG )GH ,
with appropriate si > 0. However, by increasing G the
columns of the matrix G, containing array responses over
the grid points, become more and more correlated. It is well
known from classical Compressed Sensing [26, 27] that in a
3 All the results in this paper remain valid for other grids–other than G–as far
as they are sufficiently dense. For G, the matrix G becomes an oversampled
Fourier matrix, which provides the additional advantage of reducing the
computational complexity as we explain in Section IV-D.

5

sparse estimation problem such as (7), the correlation among
the columns of the sensing matrix typically degrades the
performance of estimation of W∗ , e.g., it creates spurious
rows in W∗ . It is remarkable that, as far as estimating the
signal subspace S∗ is concerned, increasing G does not incur
any degradation of the performance, thanks to the convergence
of (7) to (5) proved in Proposition 1.
♦
In this section, we derive our low-complexity algorithm for
solving the optimization problem (7) using the well-known
Forward-Backward Splitting (FBS) for minimizing sum of two
convex functions (see [40, 41] and the refs. therein).
Definition 1: Let g : Ck → R be a convex function. The
proximal operator of g denoted by proxg : Ck → Ck is defined
♦
by proxg (x) := arg miny∈Ck g(y) + 12 kx − yk2 .
Note that for any arbitrary convex function g and a fixed x ∈
Ck , the modified convex function g(y)+ 12 kx−yk2 is strongly
convex and has a unique minimizer, thus, proxg (x) is always
well-defined (single-valued) for any arbitrary x ∈ Ck .
Consider the objective function in (7). After suitable scaling,
we can write (7) as the minimization of the convex function
f (W) = f1 (W) + f2 (W), where
1
kǦW − Xk2 , f2 (W) := kWk2,1 ,
2ζ

(8)

√
where ζ = T . The gradient of f1 is given by ∇f1 (W) =
1 H
ζ Ǧ (ǦW −X). Notice that ∇f1 is a Lipschitz function with
a Lipschitz constant β, i.e.,
k∇f1 (W) − ∇f1 (W0 )k ≤ βkW − W0 k,

1:
2:
3:
4:
5:
6:
7:

B. Forward-Backward Splitting

f1 (W) :=

Algorithm 1 FBS for l2,1 -Minimization.

(9)

with β = ζ1 λmax (ǦH Ǧ) = ζ1 λmax (ǦǦH ), where λmax
denotes the maximum singular value of a given matrix. Note
that if the grid size G is sufficiently large and the grid points
are distributed approximately uniformly over the AoAs, we
have that
G
o
1 nX
G
G
ǦǦH = B
a(θi )a(θi )H BH ≈ BIM BH = Im ,
m
m
m
i=1
G
√ . Using the standard
which implies that β = ζm
= mG
T
results, we obtain the following upper bound for f1 (W).
Proposition 2: Let W0 be a given point. Then, f1 (W) for
every W can be upper bounded by fˇ1 (W), where

(kW

k−α)

i,.
+
α > 0, (proxαf2 (W))i,. =
Wi,. is obtained by
kWi,. k
simply shrinking the rows of W, where (x)+ := max(x, 0).
With this explanation, we propose the following iterative
algorithm based on FBS. We initialize W(0) = 0 and define
for k = 1, 2, . . . the sequence W(k+1) := prox β1 f2 W(k) −

1
(k)
) . It is seen that the estimate W(k+1) at iteration
β ∇f1 (W
k is obtained by lower-bounding the function f1 (W) in
a neighborhood of W0 = W(k) by fˇ1 (W) according to
Proposition 2 and finding the optimal solution of the resulting
function fˇ(W). A variable step-size variant of our proposed
algorithm is given in Algorithm 1, where in each iteration,
the functions f1 is minimized by moving along −∇f1
with a positive step-size αk (forward step), followed by f2
being minimized by applying the proximal operator proxαk f2
(backward step). This approach is known as operator splitting
since individual components of f , i.e., f1 and f2 , are optimized
sequentially rather than jointly. The advantage is that splitting
reduces the computational complexity since most of the time
computing the joint proximal operator proxf1 +f2 is much more
complicated than computing the individual one, e.g., proxf2 .
As f (W) is strongly convex, it has a unique optimal
solution W∗ . From the convergence analysis in [40], we obtain
the following result.
Proposition 3: Let {W(k) }∞
k=0 be the sequence generated
by Algorithm 1 for an arbitrary initial point W(0) and for
arbitrary selection of step-sizes according to Algorithm 1.
∗
Then, {W(k) }∞
k=0 converges to the unique solution W . 
In order to further increase the convergence speed of
Algorithm 1, we apply Nestrov’s update rule [42], which has
been applied for the l1 -norm minimization in [43].

Algorithm 2 FBS with Nestrov’s Update.
1:
2:
3:

β
fˇ1 (W) = f1 (W0 ) + h∇f1 (W0 ), W − W0 iR + kW − W0 k2 , 4:
2
where h., .iR denotes the real-valued inner product.

Proof: Proof in Appendix B.
From Proposition 2, it follows that fˇ1 (W) gives an upper
bound on f1 (W) around a given point W0 , which is indeed
tight at W0 . This implies that f (W) can be upper-bounded
by fˇ(W) := fˇ1 (W) + f2 (W). Minimizing fˇ(W) can
be equivalently written as minimizing f2 (W) + β2 kW −
W0 + β1 ∇f1 (W0 )k2 , whose optimal solution is given by
W00 = prox β1 f2 (W0 − β1 ∇f1 (W0 )) in terms of the proximity
operator of the l2,1 -norm f2 (W) = kWk2,1 according to
Definition 1. Standard calculations show that for a given

Initialization: Fix  ∈ (0, min{1, β1 }), W(0) .
for k = 1, . . . , do
αk ∈ [, 2/β − ]
Z(k) = W(k) − αk ∇f1 (W(k) )
χk ∈ [, 1]
W(k+1) = W(k) + χk (proxαk f2 (Z(k) ) − W(k) ).
end for

5:
6:
7:
8:

Initialization: Fix W(0) , set Z(0) = W(0) , and t0 = 1.
for k = 0, 1, . . . , do
R(k) = Z(k) − β1 ∇f1 (Z(k) ).
W(k+1) = prox β1 f2 (R(k) ).
√
1+ 4t2k +1
.
tk+1 =
2
−1
αk = 1 + ttkk+1
.
Z(k+1) = W(k) + αk (W(k+1) − W(k) ).
end for

Proposition 4 (Theorem 11.3.1 in [44]): Let {W(k) }∞
k=0
be the sequence generated by Algorithm 2 for an arbitrary
initial point W(0) and for the step-sizes according to the
Nestrov’s update rule. Then, for any k, we have f (W(k+1) ) −
∗
−W(0) k2
f (W∗ ) ≤ 4βkW(k+1)
.

2

6

Remark 3: The advantage of Nestrov’s update, as seen
from Proposition 4, is that the gap to the optimal value, i.e.,
f (W(k) ) − f (W∗ ), scales like O( k12 ) as a function of the
number of iterations k rather than O( k1 ) that typically occurs
for the selection of step-sizes according to Algorithm 1. In
♦
particular, the scaling O( k12 ) is optimal [44].
Remark 4: As mentioned in Remark 2, increasing the grid
size G does not degrade the performance of the subspace
√
estimation. However, since the Lipschitz constant β = mG
T
grows proportionally to G, it is seen from Proposition 4 that
increasing G reduces the speed of the algorithm.
♦
C. Time-Varying Sampling Operators
When the dimension of sketches is much less than the
number of antennas, i.e., m  M , or when the sampling
ratio is much smaller than the normalized angular spread of
m
 2θ∆θ
, using a fixed sampling
channel vectors, i.e., M
max
matrix B results in an aliasing pattern that is not typically
resolved even by taking several measurements. Therefore, to
improve the estimation performance, it is beneficial to use
time-varying sampling matrices B(t) in each slot t ∈ [T ] such
that the aliasing caused by a sampling matrix B(t) at a time
t ∈ [T ] is resolved by other sampling matrices B(t0 ) at other
times t0 6= t. Let It ⊆ [M ] denote the indices of the sampled
antenna elements at time t ∈ [T ]. We always assume that the
indices belonging to It are sorted in an increasing ordered. We
follow the MATLAB© convention that for a vector m ∈ CM ,
we have B(t)m = m(It ), where m(It ) denotes an m-dim
vector containing the components of m belonging to It . All
the formulations for the fixed operator B can be immediately
extended to the time-variant
P case by defining f22(W) =
1
kWk2,1 and f1 (W) = 2ζ
t∈[T ] kǦt W.,t − X.,t k , where
√
ζ = T and where Ǧt = √1m B(t)G, and setting ∇f1 (W) to
be the G×T matrix with ∇f1 (W).,t = ζ1 ǦH
t (Ǧt W.,t −X.,t )
for t ∈ [T ]. Notice that ∇f1 is again a Lipschitz function in
this case with a Lipschitz constant


1
1
β = max λmax (ǦH
max λmax (Ǧt ǦH
t Ǧt ) =
t)
ζ t∈[T ]
ζ t∈[T ]
G
X

1
=
max λmax (B(t)
a(θi )a(θi )H B(t)H )
mζ t∈[T ]
i=1

G
G
G
≈
max λmax (B(t)IM B(t)H ) ≈
= √ ,
mζ t∈[T ]
ζm
m T
which is the same as in the time-invariant case. This implies
that all the steps of Algorithm 1 and 2, and their convergence
guarantee still hold in this case.
D. Computational Complexity
Each iteration of both Algorithm 1 and 2 requires computing
T columns of ∇f1 , where the t-th column, t ∈ [T ], is
given by ∇f1 (W).,t = ζ1 ǦH
t (Ǧt W.,t − X.,t ), evaluated at
W = W(k) at iteration k. For the special grid G with the
discrete AoAs θi := sin−1 (−1+ 2(i−1)
G ) sin(θmax ) , i ∈ [G],
in the angular range [−θmax , θmax ], the matrix G becomes an
oversampled Fourier matrix, namely, the columns of G are
π
c
2c
Mc T
given by (ωG
, ωG
, . . . , ωG
) , where ωG = ej G and where
c ∈ {−G, −G + 2, . . . , G − 1}. This special structure of G,

as a result that of Ǧ, can be exploited to compute ∇f1 (W)
quite efficiently.
For each t ∈ [T ], we first compute Ǧt W.,t . Following the
MATLAB© notation, let m = G ifft(W.,t , G) ∈ CG be the
inverse Discrete Fourier Transform (DFT) of W.,t scaled with
G, which can be efficiently computed using the Fast Fourier
Transform (FFT) algorithm provided that M and G are powers
of 2. Then, Ǧt W.,t is simply given by √1m m(It ), where It
denote the indices of the sampled antennas at t ∈ [T ]. The
whole complexity
 of this step of calculation for all t ∈ [T ] is
O T G log2 (G) . After computing Ǧt W.,t , t ∈ [T ], we need
to calculate ǦH
t rt , where rt = Ǧt W.,t − X.,t . This can be
simply done by setting m to be an M -dim all-zero vector,
and embedding rt in m in indices belonging to It such that
m(It ) = rt and taking the G-point DFT of m, which gives
√1 fft(m, G). The whole complexity of this step for
ǦH
t rt =
m

all t ∈ [T ] is again O T G log2 (G) .
Letting Tconv be the number of iterations necessary for
the convergence, the whole computational complexity is
O 2Tconv T G log2 (G) , which is at least two orders of
magnitude less than the complexity of directly solving the
SDP (5) with off-the-shelf SDP solvers. As we explained
in Remark 2, increasing the grid size G does not degrade
the recovery performance. However, as also mentioned in
Remark 4, it increases the Lipschitz constant β of ∇f1 and
slows down the convergence of the algorithm. The main
reason is that increasing β makes the shrinkage operation
in the proximal operator prox β1 f2 softer. As a result, the
algorithm requires more iterations to identify the dominant
grid elements. Thus, we expect that Tconv scale proportionally
G
G
to the oversampling factor M
. We always use M
= 2.
Our numerical simulations show that for this choice of
oversampling factor, both Algorithm 1 and 2, and especially
Algorithm 2, converge in only a couple of iterations.
V. E XTENSION TO OTHER A RRAY G EOMETRIES
A. 2D Rectangular Array Configurations
Our proposed algorithms can be extended to a 2D
rectangular array consisting of M = Mx My antenna elements,
arranged over a rectangular grid
R=

n

(i −

o

Mx + 1
My + 1
)dx , (j −
)dy : i ∈ [Mx ], j ∈ [My ] ,
2
2

in the 2D plane of the array, having a horizontal spacing dx and
a vertical spacing dy between its elements. We consider a 3D
Cartesian coordinate chart with an xy-plane given by the 2D
plane of the array and with a z-axis orthogonal to it. We denote
the M -dim (M = Mx My ) array responses by a(ξ), where ξ
belongs to the unit 2D sphere S2 = {ξ ∈ R3 : kξk = 1}
(lying in 3D space) and parameterizes the AoAs; sometimes,
it is better to use a coordinate chart for S2 in which every
point ξ is represented by two angles: the polar angle θ and the
azimuthal angle φ. It is more convenient to denote the M -dim
array response a(ξ) with double index (x, y) ∈ [Mx ]×[My ] :=
{(ix , iy ) : ix ∈ [Mx ], iy ∈ [My ]}, where we have
[a(ξ)]xy = ej

2π
λ hξ,rxy i

,

(10)

7


M +1
where rxy = (x − Mx2+1 )dx , (y − y2 )dy , 0 ∈ R3 denotes
the location of the array element indexed by (x, y) in the
2D plane of the array (xy-plane). The channel vector of a
user, whose scattering channel consists of a collection of
p scatterers with AoAs parameterized by {ξi : i ∈ [p]}
and
Pp channel gains {wi (t) : i ∈ [p]}, is given by h(t) =
i=1 wi (t)a(ξ
Pp i ). The channel covariance matrix is also given
by S = i=1 σi2 a(ξi )a(ξi )H , which using the double-index
notation can be represented by
[S]xy,x0 y0 =

p
X

σi2 [a(ξi )]xy [a(ξi )H ]x0 y0

ξz

ξy

ξx

Fig. 2: A non-uniform grid of AoAs over the unit sphere,
whose projection on the ξx ξy -plane is a rectangular grid.

i=1

=

p
X

σi2 ej

2π
λ hξ,rxy −rx0 y 0 i

σi2 ej

2π
λ hξ,rx−x0 ,y−y 0 i

i=1

=

p
X

.

(11)

i=1

It is seen from (11) that S has a block-Toeplitz form, i.e., it
can be represented by an M × M matrix containing Mx × Mx
blocks of dimension My × My , where the matrix at block
(x, x0 ) is given by Ux−x0 and depends only on x − x0 , where
we also have U−k = UkH , k ∈ [Mx ], due to the Hermitian
symmetry. Similarly, S can be represented with an M × M
matrix containing My × My blocks of dimension Mx × Mx ,
where the matrix at block (y, y 0 ) is given by Vy−y0 , with
V−k = VkH , k ∈ [My ]. In fact, S is even more structured
since all the diagonal blocks of S in both block representations
are equal to a Toeplitz matrix, whereas a block-Toeplitz
matrix generally might not have Toeplitz diagonal blocks.
The originally proposed AML Algorithm for the ULA in
[10, 11] can be generalized to 2D rectangular arrays. It can be
formulated as an SDP similar to (5) by replacing the constraint
set T+ with the set of PSD Hermitian block-Toeplitz matrices
denoted by BT+ , which is still a convex set.
We again assume that in each slot t ∈ [T ], we only
sample a collection of m  M array elements via a possibly
time-variant sampling matrix B(t). Similar to the previous
case for ULA, we define a 2D grid G of size G by quantizing
the continuum of AoAs, and construct the M × G matrix
consisting of the array responses over the discrete AoAs
belonging to G. A direct inspection in the proof of Proposition
1 indicates that the SDP for AML Algorithm in this case
can still be approximated by the l2,1 -norm regularized convex
optimization in (7). All the steps of the algorithm and all the
parameters remain the same as in the case of ULA. However,
due to the 2D lattice array configuration, we need to apply
2D DFT to compute ∇f1 (W) in each step rather than 1D
DFT used for the ULA. This can still be efficiently computed
provided that both Mx and My , and the oversampling ratios
Gy
Gx
Mx and My are powers of 2, where the total computational

complexity is again given by O 2Tconv T G log2 (G) .
As explained in Section IV-D, using the computational
advantage of FFT algorithm requires a special design of
the grid G that contains the AoAs θi = sin−1 (−1 +

2(i−1)
G ) sin(θmax ) , i ∈ [G]. However, if G is large enough,
G has a performance comparable with any other grid of
similar size in approximating the signal covariance matrix (see

(6)) since it covers the whole angular range [−θmax , θmax ].
Unfortunately, this is not the case for 2D rectangular arrays:
exploiting the computational advantage of 2D FFT restricts the
range of AoAs that can be processed. To explain this better,
let us consider two uniform grids: Gx a grid of size Gx in
x
x
y
y
[−ξmax
, ξmax
] and Gy a grid of size Gy in [−ξmax
, ξmax
],
y
x
y
x
where ξmax and ξmax are such that ξmax dx = ξmax dy = 1,
where we define dx = dx /(λ/2) and dy = dy /(λ/2) as the
normalized horizontal and vertical spacing between the array
x
and
elements in the 2D grid R. We also assume that ξmax
y
x
2
y
ξmax satisfy the additional constraint (ξmax ) + (ξmax )2 ≤ 1.
Let us consider the following grid consisting of G = Gx Gy
points on the unit sphere S2 , each representing a specific AoA:
G = {(ξx , ξy ,

q

1 − ξx2 − ξy2 ) : ξx ∈ Gx , ξy ∈ Gy }.

(12)

This is illustrated in Fig. 2 with the grid points lying on
the unit sphere. It is seen that G can cover only a subset
of all possible AoAs. From (12) and the definition of Gx
and Gy , it is not difficult to check that the projection of
the grid points on the ξx ξy -plane builds a 2D rectangular
x
x
y
y
grid enclosed by the rectangle [−ξmax
, ξmax
] × [−ξmax
, ξmax
].
Also, notice that the array response a(ξ) at a ξ ∈ G is
My +1
Mx +1
given by [a(ξ)]xy = ejπ(x− 2 )dx ξx ejπ(y− 2 )dy ξy , where
x ∈ [Mx ], y ∈ [My ], and ξx ∈ Gx and ξy ∈ Gy denote the
xy component of ξ. Since x, y, ξx and ξy all take values in
discrete lattices (with uniform spacing 1, 1, 1/dx and 1/dy
respectively), letting G be the M × G matrix whose columns
are given by a(ξ), ξ ∈ G, with a suitable ordering, we obtain
the 2D DFT matrix. If Mx and My are powers of 2 and
G
Gx
the oversampling ratios M
and Myy are also powers of two,
x
similarly to the case of ULA in Section IV-D, we can apply
the 2D FFT algorithm to compute ∇f1 quite fast.
x
y
Remark 5: Although ξmax
dx and ξmax
dy could be selected
to be less than 1, this unreasonably restricts the spatial
resolution of the array. For a practical design, we should
first decide on the subset of AoAs on the unit sphere that
we intend to process, with the additional constraint that the
projection of this subset on the ξx ξy -plane must lie in a
symmetric rectangular region, which is necessary in order to
take advantage of computational benefits of 2D FFT. This
x
y
yields the desired ξmax
and ξmax
. Finally, to obtain the best
spatial resolution in the desired region, the array spacings
dx and dy should be set to their maximum values such that
y
x
dx = ξmax
dy = 1.
♦
ξmax

8

B. General Array Configurations
Consider a general array configuration, in which the array
responses are parameterized by {a(ξ) : ξ ∈ Ξ}, for some
parameter set Ξ representing the AoAs. In this case, the space
of all feasible signal covariance
R matrices is given by the set of
all M × M matrices S := { Ξ γ(dξ)a(ξ)a(ξ)H : γ ∈ P(Ξ)}
where P(Ξ) denotes the space of all positive measures over
Ξ. Notice that S is indeed a convex subset (cone) of the
cone of all M × M PSD matrices. Depending on the array
geometry, it might happen that S has a simple algebraic
representation that can be exploited in the optimizations.
For example, S coincides with the space of PSD Hermitian
Toeplitz matrices T+ for the ULA, and with the space of PSD
Hermitian block-Toeplitz matrices BT+ for a 2D lattice array
configuration.
The SDP formulation (5) for the AML Algorithm can be
extended to this case by replacing T+ with S. In particular,
the equivalence between SDP (5) and the l2,1 -regularized
convex optimization (7) still holds provided that the set
Ξ is quantized with a sufficiently dense grid G, such that
supξ∈Ξ inf ξ0 ∈G ka(ξ)a(ξ)H − a(ξ 0 )a(ξ 0 )H k ≤ M holds for
a sufficiently small  ∈ (0, 1). Due to the iso-norm property
of the array response vectors {a(ξ) : ξ ∈ Ξ}, this condition
√ is
satisfied provided that supξ∈Ξ inf ξ0 ∈G ka(ξ)−a(ξ 0 )k ≤  M .
Overall, we can not expect to have the low per-iteration
computational complexity O(T G log2 (G)) obtained because
of using the FFT algorithm in the case of the ULA or the 2D
lattice array unless the covariance matrices in S have other
special algebraic structures that can be exploited to speed up
numerical computations. Furthermore, one might also need to
restrict the range of AoAs, as in the 2D lattice configuration,
to benefit this underlying algebraic structure.
VI. S UBSPACE T RACKING
A. Extending the Algorithm to the Tracking mode
Up to now, we have assumed that, although the channel
gains w(t) as in (3), and as a result the channel vectors
h(t) vary i.i.d. with time, the underlying channel geometry
{(σi2 , θP
i ) : i ∈ [p]} embedded in the covariance matrix
p
2
H
S =
remains stable for quite a long
i=1 σi a(θi )a(θi )
time, especially much longer than the window size T . This
allows the signal subspace to be estimated from the low-dim
sketches inside the window [T ], and to be used for the rest
of time. In practice, {h(t)}∞
t=1 as a stochastic process is
only locally stationary and its statistics (covariance matrix)
is piecewise constant, i.e., constant over rather long intervals
of time (time scale of one to tens of seconds) and changes with
abrupt transitions when the scattering environment of the user
changes (e.g., while moving form indoor to outdoor or turning
from one street to another for a moving vehicle). In any case,
the duration of the time intervals over which the covariance
is time-invariant is 3 up to 4 orders of magnitude larger than
the duration of the data transmission slots. Therefore, we can
collect a window of T i.i.d. samples (for a sufficiently large
T ) in the time-frequency domain over each interval [45].
Traditionally, there are two approaches in the literature to

deal with sharp transitions in signal statistics: 1) change point
detection (see [46] and refs. therein) and 2) online tracking (see
[17, 21–25] and the refs. therein). Adapted to the subspace
estimation in this paper, in the former, one applies change
point detection algorithms to identify the transition points
in the statistics, and upon identifying a transition point, the
subspace estimation algorithm is run to reestimate/update the
signal subspace from new observations. The resulting estimate
is used until the next transition point is identified. In the latter,
in contrast, upon receiving a new observation (sketch) x(t) at
time t, the tracking algorithm updates its estimate of the signal
subspace S(t) by
S(t + 1) = αS(t) + (1 − α)I (x(t)),

(13)

where α ∈ (0, 1) is the update factor, and where I (x(t))
is a subspace innovation term that depends on the newly
received sketch x(t) (see [21–25] and the refs. therein). The
choice of α makes a trade-off between the quality of the
subspace estimation in the stationary regime (variance) and
the tracking ability of the algorithm in the non-stationary
transition regime (bias)4 : the closer α to 1, the less variance
in subspace estimation in the stationary regime, and the closer
α to 0 the faster the subspace identification after occurring a
sharp transition in the non-stationary transition region. Broadly
speaking, for a given α ∈ (0, 1) the data used for subspace
estimation effectively comes from a window of Tα ≈ log(1 1 )
α
latest observations. In fact, α → 1 makes Tα larger and
improves the subspace estimation provided that the window of
observations lies in a stationary regime. However, increasing
Tα also increases the probability of having a transition in the
middle of the window, in which case the subspace estimation
algorithm requires around Tα new observations to identify the
new subspace after the transition, thus, making the tracking
algorithm less agile.
Our proposed algorithms can be run in the tracking mode
as follows. We fix a window size T , which corresponds to
selecting a suitable value for the tracking parameter α ∈ (0, 1)
in the tracking algorithm. At every time t, we always keep
the latest T sketches Wt := {x(t − T + 1), . . . , x(t)} and
update it upon receiving a new sample x(t + 1) as Wt+1 =
Wt ∪ {x(t + 1)}\{x(t − T + 1)}. We use the optimal solution
W∗ (t) of the convex optimization (7), when the matrix of
sketches is set to X = Wt , as a warm initialization to the
algorithm at time t + 1. Typically T  1, and we expect that
adding the new sketch x(t + 1) does not effect the optimal
solution considerably. In fact, for a window size T , we expect
∗
∗
(t)k2
intuitively that kW (t+1)−W
= O( T1 ). On one hand, this
kW∗ (t)k2
implies that, for a large T , only O(1) number of iterations
would be sufficient to reach from the old estimate W∗ (t) (used
as the initialization point) to the new estimate W∗ (t+1), thus,
the whole complexity of the subspace update would be of
the order O(2G log2 (G)) per each newly arrived observation.
On the other hand, this indicates that, as expected, increasing
the window size T , makes the algorithm less agile to sharp
subspace transitions since the new estimate W∗ (t+1) can not
move far from the old one W∗ (t) in a single iteration.
4 This

is well known as the bias-variance trade-off in statistics.

9

Performance of Subspace Estimator vs. SNR

B. Further Simplified Subspace Tracking

VII. S IMULATION R ESULTS
In this section, we evaluate the performance of our
proposed subspace estimation/tracking algorithm empirically
via numerical simulations.
A. Array Model
We consider a ULA with M = 64 antennas, where in each
training period, we randomly sample only m = 16 of them,
m
= 0.25. We assume that
thus, a sampling ratio of ρ := M
the array has θmax = 60 degrees and scans an angular range
of ∆θmax = 2θmax = 120 degrees. For all the simulations,
we use a grid size of G = 2M , thus, a grid oversampling
factor of 2, and optimize (7) by running Algorithm 2 with the
Nestrov’s update rule and with i.i.d. time-varying sampling
matrices explained in Section IV-C.
B. Scaling with respect to Training Signal-to-Noise Ratio
We consider a scattering geometry, in which the received
signal power of a given user is uniformly distributed over the
angular range Θ = [10, 30] degrees, with an angular spread
of 20 degrees. We define the training Signal-to-Noise Ratio
2
]
(SNR) by E[kh(t)k
E[kn(t)k2 ] , where h(t) denotes the user channel
vector and where n(t) is the array noise at time t as in (3).
Fig. 3 illustrates the scaling of the performance metric
b ) defined in Section II-C versus the training SNR for
Γ(p, p
different training lengths (window sizes) T ∈ {50, 100, 200}.
To obtain the curve for each T , we average the resulting
performance of the subspace estimator versus SNR over 100
independent simulations, where in each simulation we run the
algorithm for approximately 50 iterations before convergence.
We consider two different family of measurement matrices for
simulation: i) The 0-1 sampling matrices denoted by “Bin”,
where each row of the matrix contains one 1 at a specific
random column and corresponds to a random antenna selection
as explained in Section II-B, and ii) The random phase-shift
matrices denoted by “PS”, where the sampling matrices B(t)
jθr,c (t)
are generated according to [B(t)]r,c = e √M , where θr,c (t)
are selected i.i.d. randomly in each row r, column c, and
across different slots t ∈ [T ] form the set of quantized phases

1
b)
Performance Metric Γ(p, p

In practical implementations, updating the subspace at each
time t requires the following steps: updating the weighting
matrix W(t), computing the estimate of signal covariance
matrix according to Proposition 1, and computing the SVD of
the resulting covariance matrix and identifying its dominant
subspace. This might be too complicated in some real-time
implementations in massive MIMO. Instead, we can use the
weighting matrix W(t) to identify the position of dominant
elements (active elements) in the over-complete dictionary
over the grid given by G. This only requires updating the
l2 -norm or equivalently the l2 -norm squared of the rows of
W(t), i.e., kW(t)i,. k2 for i ∈ [G], after each iteration t, which
can be done quite fast. At each time t, the M × q submatrix of
G (for some q  M ) corresponding to the dominant rows of
W(t) with significantly large l2 norms provides an estimate
of the dominant signal subspace.

0.8
0.6

Bin
PS
Bin
PS
Bin
PS

0.4
0.2

0
−10

−5

0

5

T
T
T
T
T
T

= 50
= 50
= 100
= 100
= 200
= 200

10

SNR [dB]

Fig. 3: Performance of our proposed subspace estimation algorithm
versus SNR for different training lengths T ∈ {50, 100, 200} for 0-1
sampling matrices “Bin” and random phase-shift matrices “PS”. The
b ) 9 1 for large SNR is due to the finite grid size G.
reason Γ(p, p

: k = 0, . . . , 2b − 1} with b = 5 bit precision. In
{ k2π
2b
massive MIMO applications, the former family of sampling
matrices can be implemented via switches whereas the latter
is more suitable for implementation via constant-amplitude
phase-shafting networks. It is seen from Fig. 3 that both family
of sampling matrices have quite similar performances although
the 0-1 sampling matrices yield a more efficient numerical
implementation via the FFT algorithm as explained in Section
IV-D by essentially avoiding any matrix multiplication. The
simulation results are qualitatively similar to the results in
[10, 11] but are obtained for a different performance metric.
In applications such as JSDM, the practically important SNR
regime is around 0 dB up to 10 dB in which the system has
a considerably high throughput (measured in terms of the
achievable sum-rate or spectral efficiency). Fig. 3 illustrates
that in this regime of SNR, our proposed subspace estimation
scheme has an excellent performance.
C. Comparison with SVT Algorithm
We compare the performance of our algorithm with that
of the SVT (singular value thresholding) Algorithm in [38],
which provides the state of the art performance in subspace
estimation (and matrix completion). SVT is a batch algorithm,
i.e., it estimates the signal subspace from a fixed number
of sketches T , in contrast with an online tracking algorithm
where the number of sketches increases with time5 . Each
iteration of SVT consists of computing the SVD of an M × T
matrix followed by a singular value thresholding, which can
be done efficiently (e.g., via Lanczos algorithm) when the
matrix is quite low-rank. However, our algorithm is much
faster since it does not requires any SVD computation or
matrix multiplication as all the calculations are done via the
FFT algorithm.
For simulations, we assume that the received signal power
5 Although the effective number of sketches still remains constant and
depends on the subspace update factor α as in (13).

10

Performance of Subspace Estimator vs. SNR

b)
Performance Metric Γ(p, p

1
b)
Performance Metric Γ(p, p

SNR= −5 dB

1

0.8

0.8

0.8

0.6

0.6

0.4
0.2

0.6

0

0

0.4

T = 50
T = 100
T = 200
100

T = 50
T = 100
T = 200

0.2
200

300

400

0

0

100

200

Time

0.4

0

5

10

T
T
T
T

15

= 100,
= 400,
= 800,
= 1600,

20

ρ = 0.25
ρ = 0.5
ρ = 0.5
ρ = 0.5

25

SNR= 5 dB

30

SNR [dB]

Fig. 4: Comparison of performance versus SNR of our algorithm
(Alg. 2) for with that of SVT for different number of sketches T and
m
.
different sampling ratios ρ = M

of the user is uniformly distributed over the angular range
Θ = [10, 30] degrees as in Section VII-B. Fig. 4 illustrates the
comparison of the performance versus SNR of our algorithm
for T = 100 and a sampling ratio ρ = 0.25 with that of
SVT for different T ∈ {400, 800, 1600} and a sampling ratio
ρ = 0.5. It is seen that our algorithm performs much better
than SVT even for a smaller data size T and antenna sampling
ratio ρ. In particular, for a target estimation performance, it
runs orders of magnitude faster since it requires smaller T and
consumes much less storage (scaled proportionally to M ×T ).
D. Tracking Performance
Fig. 5 illustrates the tracking performance of our proposed
estimator for different window size T ∈ {50, 100, 200} and
different SNR when there is a sharp transition in the channel
statistics (geometry). We consider a time interval of length
400 for simulation that contains a sharp transition in the
middle at time ttr = 200. For t = 1, . . . , ttr − 1, the user
signal power is uniformly distributed in the angular range
Θ = [10, 30] degrees, with an angular spread of 20 degrees
as in Section VII-B, whereas at time t = ttr , the geometry
of the channel changes abruptly such that for t = ttr , . . . , 400
the user signal power is uniformly distributed in the angular
range Θ0 = [−40, −20], where Θ ∩ Θ0 = ∅.
Fig. 5 illustrates a random sample path of the performance
b ) generated by the algorithm for the 0-1 sampling
metric Γ(p, p
matrices during the whole simulation. To generate these plots,
we run our proposed subspace estimation algorithm in a
tracking mode as explained in Section VI, in which upon
receiving a new sketch, we run only one iteration of our
proposed algorithm while treating the previous estimate as
initialization. We start the algorithm with the zero initialization
at time t = 0, where it is seen from Fig. 5 that the algorithm
identifies the signal subspace in quite a short time. It is also
seen that immediately after the sharp transition in the channel,
b ) → 0, however, the algorithm
the performance metric Γ(p, p
is able to track/identify the new signal subspace in quite a short

0.8

0.8

0.6

0.6

0.4

0

0

400

300

400

SNR= 10 dB
1

0.2

300

Time

1

b)
Performance Metric Γ(p, p

Alg. 2
SVT
SVT
SVT

0.2

0
−5

SNR= 0 dB
1

0.4

T = 50
T = 100
T = 200
100

T = 50
T = 100
T = 200

0.2
200

300

400

0

0

100

Time

200

Time

Fig. 5: Tracking performance of our proposed algorithm (Alg. 2) for
different window sizes T ∈ {50, 100, 200} and for different training
SNR. There is a transition in statistics at time ttr = 200.

time. Interestingly, it is seen that for a window of size T , the
delay before identifying/tracking the new subspace is around
T
2 , namely, immediately after half the observation window
is filled with new sketches generated with the new channel
geometry, the algorithm makes a sharp transition from the old
signal subspace to the new one.
E. Comparison with PETRELS Algorithm
We compare the performance of our algorithm in the
tracking mode with that of PETRELS Algorithm proposed in
[25]. Note that PETRELS is an online algorithm (in contrast
with SVT which is a batch algorithm) and provides the state
of the art performance in online subspace estimation. In the
online mode the number of sketches is not fixed (in contrast
with the batch mode) and increases with time. Denoting by
T the size of the window of sketches used by our proposed
algorithm in the tracking mode, we set the update factor of
κ
κ
PETRELS to e− T , i.e., α = e− T with our notation in (13),
where we tune κ ∈ [1, 5] to obtain the best performance. In
this way, we make sure that the effective number of sketches
used by PETRELS comes from a window of approximate size
T of the latest sketches.
Fig. 6 illustrates the simulation results. We compare the
performance of PETRELS with that of our proposed tracking
algorithm for different T and antenna sampling ratio ρ. It
is evidently seen that our algorithm has a much superior
performance in the stationary regime and tracks the subspace
transitions much faster in the non-stationary regime. Fig. 6 also
illustrates the interesting bias-variance trade-off underlying the
PETRELS performance (which also exists in our proposed
tracking algorithm as can be seen from Fig. 5), where for
a fixed sampling ratio ρ, the performance of PETRELS
in the stationary regime (variance) improves by increasing
the window size T , but this comes at the price of much

11

Tracking active elements for SNR=10 dB and T = 50
−1

Performance of Subspace Tracking for SNR= 10 dB

b)
Performance Metric Γ(p, p

1

sin(θ)/sin(θmax )

0.8
0.6
0.4

Alg. 2
PET.
PET.
PET.
PET.

0.2
0

0

0.2

0.4

T
T
T
T
T

= 100, ρ = 0.25
= 400, ρ = 1
= 800, ρ = 1
= 400, ρ = 0.66
= 800, ρ = 0.66

0.6
Time

0.8

−0.5
0
0.5
1

1
·104

Fig. 6: Comparison of tracking performance of our algorithm (Alg. 2)
with that of PETRELS for SNR= 10 dB.

slower tracking in the non-stationary regime. It is also seen
from Fig. 6 that the performance of PETRELS degrades
considerably by reducing the antenna sampling factor ρ,
whereas our algorithm is less sensitive to ρ and works perfectly
even with a quite low ρ = 0.25.
F. Simplified Subspace Tracking
As explained in Section VI-B, a more low-complexity
estimate of the signal subspace at each time t can be obtained
by identifying the dominant rows of the weighting matrix
W(t). Fig. 7 illustrates the strength of different rows of W(t)
at time t, corresponding to the estimated received power from
different angular grid element in a tracking mode during
t = 1, . . . , 400. We assume that as in the simulations in
Section VII-D, the received angular power distribution of
the user undergoes a sharp transition from the angular range
Θ = [10, 30] to Θ0 = [−40, −20] at time ttr = 200. It is seen
from Fig. 7 that although there are some spurious rows, our
proposed algorithm tracks the location (support) of dominant
rows corresponding to the AoAs of the scatterers very well.
VIII. C ONCLUSION
In this paper, we proposed an efficient and low-complexity
subspace estimation algorithm, with a special focus towards
massive MIMO applications. We mainly studied the AML
Algorithm proposed in [10, 11], where we showed that the
quite slow and time-consuming SDP optimization of AML
Algorithm in [10, 11] (especially when the antenna size
M is quite large) can be well approximated with another
convex optimization problem, for which we derived a novel
iterative low-complexity algorithm. We also considered a
generalization of the original AML Algorithm in which the
projection (sampling) operator may be time-variant, resulting
in a further improvement in subspace estimation. We explained
how our proposed algorithm can be applied to more practical
array configurations such as 2D rectangular lattice arrays and
provided guidelines for efficient numerical implementation for
general array configurations. We also extended our proposed

0

100

200
Time

300

400

Fig. 7: An image of the strength of different grid elements
estimated by our proposed algorithm during a tracking period
t = 1, . . . , 400. The angular power distribution of the user
makes a transition from Θ = [10, 30] to Θ0 = [−40, −20] at
time ttr = 200.
algorithm such that it can be run in the online tracking mode.
We evaluated the estimation/tracking performance of our
algorithm empirically via numerical simulations and compared
it with the performance of other state of the art subspace
estimation/tracking algorithms in the literature.
A PPENDIX A
P ROOF OF P ROPOSITION 1
The proof follows by extending Theorem 1 in [47]. The key
observation is that for a vector w ∈ CT , the l2 -norm kwk of
w can be written as the output of the following optimization
kwk =

kvk2 + |s|2
.
2
v∈CT , s∈C: sv=w
min

(14)

In particular, denoting by (v∗ , s∗ ) the optimal solution of (14),
we have that kwk = |s∗ |2 . Applying this to the rows of a G×T
coefficient matrix W as in (7), we obtain that
kWk2,1 =

kVk2 + kΓk2
,
2
V∈CG×T , Γ∈D: ΓV=W
min

(15)

where D denotes the space of G × G diagonal matrices with
diagonal elements in C, and where Γ = diag(γ1 , . . . , γG ) ∈ D.
In particular, similar to (14) we have that kWi,. k = |γi∗ |2 ,
∗
where Γ∗ = diag(γ1∗ , . . . , γG
) is the optimal solution of (15).
Replacing kWk2,1 in (7) with (15), we can transform (7) into
the following optimization problem
(V∗ , Γ∗ ) =

arg min
V∈CG×T , Γ∈D

kǦΓV − Xk2
√
+ kVk2 + kΓk2 .
T

For a fixed Γ, the minimizing matrix V as a function of Γ
can be obtained via a least-square minimization, where after
replacing the solution and applying the matrix inversion lemma
[48] and further simplifications, we obtain the following
optimization in terms of Γ
H
 ΓΓH

b x + tr( ΓΓ
√ ).
Γ∗ = arg min tr (Ǧ √ ǦH + Im )−1 C
T
T
Γ∈D

12

This optimization can be reparameterized with P =
2

H
ΓΓ
√
T

=

diag(p1 , . . . , pG ) ∈ D+ , where pi = |γ√iT| ∈ R+ , for i ∈ [G],
and where D+ denotes the set of all G × G diagonal matrices
with positive diagonal elements. With this parametrization, we
obtain the following optimization for the matrix P ∈ D+ :


b x + tr(P),
P∗ = arg min tr (ǦPǦH + Im )−1 C
(16)
P∈D+

Moreover, denoting by P∗ = diag(p∗1 , . . . , p∗G ) the solution of
(16), we have the following relation
∗
kWi,.
k
|γ ∗ |2
,
p∗i = √i = √
T
T

(17)

1
1
GPGH = G diag(p1 , . . . , pG )GH ,
(18)
m
m
for some appropriate P ∈ D+ with pi ≥ 0, i ∈ [G]. This
implies that the term ǦPǦH in (16) can be replaced with
BSBH (recall that Ǧ = √1m BG), where S takes values from
the convex set of all feasible signal covariance matrices, which
in the case of a ULA coincides with T+ . Also note that due
to 0-1 sampling matrix B and the special structure of the
array responses, every column of Ǧ has a unit l2 -norm, which
implies that
S≈

G
X
i=1

pi kǧi k2 =

G
X

pi = tr(P),

(19)

i=1

where ǧi denotes the i-th column of Ǧ, which has a unit norm.
Again replacing ǦPǦH by BSBH , it results that (16) is well
approximated by the following convex optimization


b x + tr(BSBH ). (20)
S∗ = arg min tr (BSBH + Im )−1 C
S∈T+

Using the well-known Schur’s complement condition for
positive semi-definiteness (see [49] page 28), we can write (20)
in the form of SDP (5) for the AML Algorithm as in [10, 11].
In particular, having the optimal solution W∗ of (7), or the
optimal solution P∗ of (16), from the parametrization in (18),
the optimal solution S∗ of the SDP (5) can be approximated
kW∗ k
p∗
by S∗ = G diag(s∗1 , . . . , s∗G )GH , where s∗i = mi = m√i,.T ,
for i ∈ [G], and where we used (17). This completes the proof.
A PPENDIX B
P ROOF OF P ROPOSITION 2

Let s ∈ [0, 1] and let us define ∆(s) := W0 + s(W − W0 )
and h(s) := f1 (∆(s)). We have
f1 (W) − f1 (W0 ) = f1 (∆(1)) − f1 (∆(0))
Z 1
= h(1) − h(0) =
h0 (s)ds
0
Z 1
=
h∇f1 (∆(s)), W − W0 iR ds
0

0

(a)

≤ h∇f1 (W0 ), W − W0 iR
Z 1
+
βk∆(s) − W0 kkW − W0 kds
0
Z 1
0
0
= h∇f1 (W ), W − W iR +
sβkW − W0 k2 ds
0

β
= h∇f1 (W ), W − W iR + kW − W0 k2 ,
2
0

between the optimal solution W∗ of (7) and the optimal
solution P∗ = diag(p∗1 , . . . , p∗G ) of (16). Note that as in
Section IV-A, we assume that the grid G is dense enough such
that any signal covariance matrix can be well approximated by

tr(ǦPǦH ) =

= h∇f1 (W0 ), W − W0 iR
Z 1
+
h∇f1 (∆(s)) − ∇f1 (W0 ), W − W0 iR ds

0

where in (a) we used the Cauchy-Schwarz inequality and the
Lipschitz property of ∇f1 . This completes the proof.
R EFERENCES
[1] T. L. Marzetta, “Noncooperative cellular wireless with unlimited
numbers of base station antennas,” IEEE Trans. on Wireless
Commun., vol. 9, no. 11, pp. 3590–3600, Nov. 2010.
[2] H. Huh, G. Caire, H. Papadopoulos, and S. Ramprashad,
“Achieving massive MIMO spectral efficiency with a
not-so-large number of antennas,” IEEE Trans. on Wireless
Commun., vol. 11, no. 9, pp. 3226–3239, 2012.
[3] 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. on Sel. Areas on Commun. (JSAC), vol. 31, no. 2, pp.
160–171, 2013.
[4] E. Larsson, O. Edfors, F. Tufvesson, and T. Marzetta,
“Massive MIMO for next generation wireless systems,” IEEE
Communications Magazine, vol. 52, no. 2, pp. 186–195, 2014.
[5] C. Shepard, H. Yu, N. Anand, E. Li, T. Marzetta, R. Yang,
and L. Zhong, “Argos: Practical many-antenna base stations,”
in Proceedings of the 18th Annual International Conference on
Mobile Computing and Networking. ACM, 2012, pp. 53–64.
[6] R. Rogalin, O. Y. Bursalioglu, H. Papadopoulos, G. Caire, A. F.
Molisch, A. Michaloliakos, V. Balan, and K. Psounis, “Scalable
synchronization and reciprocity calibration for distributed
multiuser MIMO,” IEEE Trans. on Wireless Commun., vol. 13,
no. 4, pp. 1815–1831, 2014.
[7] A. Adhikary, J. Nam, J.-Y. Ahn, and G. Caire, “Joint spatial
division and multiplexing: the large-scale array regime,” IEEE
Trans. on Inform. Theory, vol. 59, no. 10, pp. 6441–6463, 2013.
[8] J. Nam, A. Adhikary, J.-Y. Ahn, and G. Caire, “Joint spatial
division and multiplexing: Opportunistic beamforming, user
grouping and simplified downlink scheduling,” IEEE J. of Sel.
Topics in Sig. Proc. (JSTSP), vol. 8, no. 5, pp. 876–890, 2014.
[9] A. Adhikary, E. Al Safadi, M. K. Samimi, R. Wang, G. Caire,
T. S. Rappaport, and A. F. Molisch, “Joint spatial division and
multiplexing for mm-wave channels,” IEEE J. on Sel. Areas on
Commun. (JSAC), vol. 32, no. 6, pp. 1239–1255, 2014.
[10] S. Haghighatshoar and G. Caire, “Channel vector subspace
estimation from sample covariance of low-dimensional
projections,”
in
International
Zurich
Seminar
on
Communications, 2016, p. 103.
[11] ——, “Massive MIMO channel subspace estimation from
low-dimensional projections,” IEEE Transactions on Signal
Processing, vol. 65, pp. 303–318, 2017.
[12] R. Kumaresan and D. W. Tufts, “Estimating the angles of arrival
of multiple plane waves,” IEEE Transactions on Aerospace and
Electronic Systems, no. 1, pp. 134–139, 1983.
[13] R. Roy and T. Kailath, “Esprit-estimation of signal parameters
via rotational invariance techniques,” IEEE Transactions on
Acoustics, Speech and Signal Processing, vol. 37, no. 7, pp.
984–995, 1989.
[14] R. O. Schmidt, “Multiple emitter location and signal parameter

13

[15]
[16]

[17]
[18]
[19]
[20]
[21]
[22]
[23]
[24]

[25]

[26]
[27]

[28]
[29]
[30]

[31]
[32]

estimation,” IEEE Transactions on Antennas and Propagation,
vol. 34, no. 3, pp. 276–280, 1986.
S. Shahbazpanahi, S. Valaee, and M. H. Bastani, “Distributed
source localization using esprit algorithm,” IEEE Transactions
on Signal Processing, vol. 49, no. 10, pp. 2169–2178, 2001.
D. W. Stein, S. G. Beaven, L. E. Hoff, E. M. Winter,
A. P. Schaum, and A. D. Stocker, “Anomaly detection from
hyperspectral imagery,” IEEE signal processing magazine,
vol. 19, no. 1, pp. 58–69, 2002.
A. H. Sayed, Fundamentals of adaptive filtering. John Wiley
& Sons, 2003.
D. Tse and P. Viswanath, Fundamentals of wireless
communication. Cambridge university press, 2005.
I. T. Jolliffe, “Principal component analysis and factor analysis,”
in Principal component analysis. Springer, 1986, pp. 115–128.
A. Papoulis and S. U. Pillai, Probability, random variables, and
stochastic processes. Tata McGraw-Hill Education, 2002.
P. Comon and G. H. Golub, “Tracking a few extreme singular
values and vectors in signal processing,” Proceedings of the
IEEE, vol. 78, no. 8, pp. 1327–1343, 1990.
B. Yang, “Projection approximation subspace tracking,” IEEE
Transactions on Signal Processing, vol. 43, no. 1, pp. 95–107,
1995.
K. Crammer, “Online tracking of linear subspaces,” Lecture
notes in computer science, vol. 4005, p. 438, 2006.
L. Balzano, R. Nowak, and B. Recht, “Online identification and
tracking of subspaces from highly incomplete information,” in
48th Annual Allerton Conference on Communication, Control,
and Computing (Allerton), 2010. IEEE, 2010, pp. 704–711.
Y. Chi, Y. C. Eldar, and R. Calderbank, “Petrels: Parallel
subspace estimation and tracking by recursive least squares from
partial observations,” IEEE Transactions on Signal Processing,
vol. 61, no. 23, pp. 5947–5959, 2013.
D. L. Donoho, “Compressed sensing,” IEEE Transactions on
Information Theory, vol. 52, no. 4, pp. 1289–1306, 2006.
E. J. Candes and T. Tao, “Near-optimal signal recovery
from random projections: Universal encoding strategies?” IEEE
Transactions on Information Theory, vol. 52, no. 12, pp.
5406–5425, 2006.
E. J. Candès and B. Recht, “Exact matrix completion via convex
optimization,” Foundations of Computational mathematics,
vol. 9, no. 6, pp. 717–772, 2009.
E. J. Candès and T. Tao, “The power of convex relaxation:
Near-optimal matrix completion,” IEEE Transactions on
Information Theory, vol. 56, no. 5, pp. 2053–2080, 2010.
E. J. Candès and C. Fernandez-Granda, “Towards a
mathematical theory of super-resolution,” Communications on
Pure and Applied Mathematics, vol. 67, no. 6, pp. 906–956,
2014.
R. Baraniuk and P. Steeghs, “Compressive radar imaging,” in
IEEE Radar Conference, 2007. IEEE, 2007, pp. 128–133.
M. Herman, T. Strohmer et al., “High-resolution radar via
compressed sensing,” IEEE Transactions on Signal Processing,

vol. 57, no. 6, pp. 2275–2284, 2009.
[33] W. U. Bajwa, J. Haupt, A. M. Sayeed, and R. Nowak,
“Compressed channel sensing: A new approach to estimating
sparse multipath channels,” Proceedings of the IEEE, vol. 98,
no. 6, pp. 1058–1076, 2010.
[34] J. A. Tropp, A. C. Gilbert, and M. J. Strauss, “Algorithms
for simultaneous sparse approximation. part i: Greedy pursuit,”
Signal Processing, vol. 86, no. 3, pp. 572–588, 2006.
[35] K. Lee, Y. Bresler, and M. Junge, “Subspace methods for joint
sparse recovery,” IEEE Transactions on Information Theory,
vol. 58, no. 6, pp. 3613–3641, 2012.
[36] J. A. Tropp, “Algorithms for simultaneous sparse
approximation. part ii: Convex relaxation,” Signal Processing,
vol. 86, no. 3, pp. 589–602, 2006.
[37] D. Malioutov, M. Çetin, and A. S. Willsky, “A sparse signal
reconstruction perspective for source localization with sensor
arrays,” IEEE Transactions on Signal Processing, vol. 53, no. 8,
pp. 3010–3022, 2005.
[38] J.-F. Cai, E. J. Candès, and Z. Shen, “A singular value
thresholding algorithm for matrix completion,” SIAM Journal
on Optimization, vol. 20, no. 4, pp. 1956–1982, 2010.
[39] A. W. Marshall, I. Olkin, and B. C. Arnold, Inequalities: theory
of majorization and its applications. Springer, 1979, vol. 143.
[40] P. L. Combettes and V. R. Wajs, “Signal recovery by
proximal forward-backward splitting,” Multiscale Modeling &
Simulation, vol. 4, no. 4, pp. 1168–1200, 2005.
[41] P. L. Combettes and J.-C. Pesquet, “Proximal splitting methods
in signal processing,” in Fixed-point algorithms for inverse
problems in science and engineering.
Springer, 2011, pp.
185–212.
[42] Y. Nesterov, “A method of solving a convex programming
problem with convergence rate O(1/k2 ),” in Soviet
Mathematics Doklady, vol. 27, no. 2, 1983, pp. 372–376.
[43] A.
Beck
and
M.
Teboulle,
“A
fast
iterative
shrinkage-thresholding algorithm for linear inverse problems,”
SIAM journal on imaging sciences, vol. 2, no. 1, pp. 183–202,
2009.
[44] A. Nemirovski, “Efficient methods in convex programming,”
2005.
[45] K. Mahler, W. Keusgen, F. Tufvesson, T. Zemen, and G. Caire,
“Propagation of multipath components at an urban intersection,”
in Vehicular Technology Conference (VTC Fall), 2015 IEEE
82nd. IEEE, 2015, pp. 1–5.
[46] H. V. Poor and O. Hadjiliadis, Quickest detection. Cambridge
University Press Cambridge, 2009, vol. 40.
[47] C. Steffens, M. Pesavento, and M. E. Pfetsch, “A compact
formulation for the `2,1 mixed-norm minimization problem,”
arXiv preprint arXiv:1606.07231, 2016.
[48] W. W. Hager, “Updating the inverse of a matrix,” SIAM review,
vol. 31, no. 2, pp. 221–239, 1989.
[49] S. P. Boyd, L. El Ghaoui, E. Feron, and V. Balakrishnan, Linear
matrix inequalities in system and control theory. SIAM, 1994,
vol. 15.

