IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. XX, NO. XX, MONTH, 2018

1

Superfast Line Spectral Estimation

arXiv:1705.06073v2 [eess.SP] 15 Feb 2018

Thomas L. Hansen, Bernard H. Fleury, Senior Member, IEEE, and Bhaskar D. Rao, Fellow, IEEE

Abstract—A number of recent works have proposed to solve
the line spectral estimation problem by applying off-the-grid
extensions of sparse estimation techniques. These methods are
preferable over classical line spectral estimation algorithms
because they inherently estimate the model order. However, they
all have computation times which grow at least cubically in
the problem size, thus limiting their practical applicability in
cases with large dimensions. To alleviate this issue, we propose
a low-complexity method for line spectral estimation, which also
draws on ideas from sparse estimation. Our method is based on
a Bayesian view of the problem. The signal covariance matrix
is shown to have Toeplitz structure, allowing superfast Toeplitz
inversion to be used. We demonstrate that our method achieves
estimation accuracy at least as good as current methods and that
it does so while being orders of magnitudes faster.

I. I NTRODUCTION
The problem of line spectral estimation (LSE) has received
significant attention in the research community for at least 40
years. The reason is that many fundamental problems in signal
processing can be recast as LSE; examples include direction
of arrival estimation using sensor arrays [1], [2], bearing
and range estimation in synthetic aperture radar [3], channel
estimation in wireless communications [4] and simulation of
atomic systems in molecular dynamics [5].
In trying to solve the LSE problem, classical approaches
include subspace methods [6] such as MUSIC [7] or ESPRIT
[8] which estimate the frequencies based on an estimate
of the signal covariance matrix. These approaches must be
augmented with a method for estimation of the model order.
Popular choices include generic information theoretic criteria
(e.g. AIC, BIC) or more specialized methods, such as SORTE
[9] which is based on the eigenvalues of the estimated signal
covariance matrix. Subspace methods typically perform extremely well if the model order is known, but their estimation
accuracy can degrade significantly if the model order is
unknown.
The stochastic maximum likelihood (ML) method is known
to be asymptotically efficient (it attains the Cramér-Rao bound
as the problem size tends to infinity) [2]. Unfortunately it also
requires knowledge of the model order.
Inspired by the ideas of sparse estimation and compressed
sensing, many papers on sparsity-based LSE algorithms have
T. L. Hansen and B. H. Fleury are with the Department of Electronic
Systems, Aalborg University, Aalborg, Denmark. B. D. Rao is with the
Electrical and Computer Engineering department, University of California,
San Diego.
The work of T. L. Hansen is supported by the Danish Council for
Independent Research under grant id DFF–4005–00549.
c 2018 IEEE. Personal use of this material is permitted. Permission from
IEEE must be obtained for all other uses, in any current or future media,
including reprinting/republishing this material for advertising or promotional
purposes, creating new collective works, for resale or redistribution to servers
or lists, or reuse of any copyrighted component of this work in other works.
Digital Object Identifier 10.1109/TSP.2018.2807417

appeared in recent years, e.g. [1], [10]. In particular, the
LSE problem is simplified to a finite sparse reconstruction
problem by restricting the frequencies to a grid. Such methods
inherently estimate the model order, alleviating the issues
arising from separate model order and frequency estimation
in classical methods. The granularity of the grid leads to
a non-trivial trade-off between accuracy and computational
requirements. To forego the use of a grid, so-called off-the-grid
compressed sensing methods have been proposed [11]–[13].
These methods provably recover the frequencies in the noisefree case under a minimum separation condition. They suffer
from prohibitively high computational requirements even for
moderate problem sizes, see Sec VI.
In [14]–[16] a Bayesian view is taken on the LSE problem. The model used in stochastic ML is extended with a
sparsity-promoting prior on the coefficients of the sinusoid
components. Thereby inherent estimation of the model order
is achieved. These algorithms generally have high estimation
accuracy. Their per-iteration computational complexity is cubic in the number of sinusoidal components, meaning that their
runtime grows rapidly as the number of components increases.
In this work we introduce the Superfast LSE algorithm for
solving the LSE problem in scenarios where the full measurement vector is available (complete data case). The modelling
and design of the basic algorithm which we present in Sec. II is
based upon the ideas in [14]–[16]. The main novelty resides
in the computational aspects of Superfast LSE. The derived
method is based upon several techniques: a so-called superfast
Toeplitz inversion algorithm [17], [18] (thereof the name of
our algorithm), low-complexity Capon beamforming [19], the
Gohberg-Semencul formula [20] and non-uniform fast Fourier
transforms [21], [22]. The Superfast LSE algorithm has the
following virtues: It inherently estimates all model parameters
such as the noise variance and model order and it has low
per-iteration computational complexity. Specifically it scales
as O(N log2 N ) where N is the length of the observed vector.
We show empirically that it converges after a few iterations
(typically less than 20). This means that for large problem sizes
our algorithm can have computation time orders of magnitude
lower than that of current methods. It does so without any
penalty in estimation accuracy. Our numerical experiments
show that Superfast LSE has high estimation accuracy across
a wide range of scenarios, being on par with or better than
state-of-the-art algorithms.
Synergistically and computationally efficiently combining
the steps in the algorithm might appear easy after the fact.
This is however not the case. Some other LSE algorithms can
benefit in terms of computational effort from our approach, yet
not to the extent achieved with the proposed algorithm. For
instance, the computational methods in Sec. III can be embedded in VALSE [16]. The resulting scheme will have high

2

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. XX, NO. XX, MONTH, 2018

computational complexity due to the variational estimation
of the posterior on the frequencies. Note that our algorithm
performs on par with VALSE, but at a significantly reduced
computational effort.
For completeness we also present a semifast version of the
algorithm which works when only a subset of entries in the
measurement vector are available. The Semifast LSE algorithm
has per-iteration complexity O(N K̂ 2 + N log N ), where K̂ is
the number of estimated sinusoids. Algorithms with similar
per-iteration complexity are derived in [14], [15], [23], [24].
We have observed that our algorithm converges in a smaller
number of iterations when compared to the algorithm in [15],
thus leading to lower total runtime.
Outline: In Sec. II we present our modelling and algorithm
for LSE. Our low-complexity computational methods are
presented in Sec. III (complete date case) and IV (incomplete
data case). In Sec. V the algorithm is extended to the case
of multiple measurement vectors. Numerical experiments are
presented in Sec. VI and conclusions are given in Sec. VII.
Notation: We write vectors as a and matrices as A. The
ith entry of vector a is denoted ai or [a]i ; the i, jth entry of
matrix A is denoted Ai,j . Let b be a binary vector (containing
only zeros and ones) of the same dimension as a, then ab
denotes a vector which contains those entries in a where the
corresponding entry in b is one. The Hadamard (entrywise)
product is denoted by .

estimation of the model order K, we follow [14], [16] and employ a model with Kmax ≥ K components1 . Each component
has an associated activation variable zk ∈ {0, 1} which is set
to 0 or 1 to deactivate or activate it. The activation variables
are collected in the sparse vector z. The effective estimated
model order is given by the number of active components.
Based on (1) we write our estimation model
y=

K
max
X

Φψ(θk )zk αk + w = A(θz )αz + w,

(2)

k=1

where θk ∈ [0, 1) and αk ∈ C are frequencies and coefficients
for k = 1, . . . , Kmax and we have defined A(θ) , ΦΨ(θ).
Due to the Gaussian noise assumption we have
p(y|α, z, θ; β) = CN(y; A(θz )αz , βI),

(3)

where CN(y; µ, Σ) denotes the probability density function
of a circularly symmetric complex normal random variable
y with mean µ and covariance matrix Σ. We assume β ∈
[εβ , ∞), where εβ > 0 is an arbitrarily small constant which
guarantees that the likelihood function is bounded below. A
Bernoulli prior is used to promote deactivation of some of the
components:
p(z; ζ) =

K
max
Y

ζ zk (1 − ζ)1−zk ,

(4)

k=1

II. A N A LGORITHM FOR L INE S PECTRAL E STIMATION
We now detail the observation model and the specific
objective of the LSE problem. The observation vector y ∈ CM
contains time-domain samples and is given by
y=

K
X

where ζ ∈ [0, 1/2] is the activation probability. The restriction
ζ ≤ 1/2 ensures that the prior is sparsity inducing. The coefficients are assumed to be independent zero-mean Gaussian
p(α; γ) =

Φψ(θ̃k )α̃k + w = ΦΨ(θ̃)α̃ + w,

where the steering vector function ψ(θk ) : [0, 1) → CN ×1
gives a Fourier vector, i.e., it has nth entry [ψ(θk )]n ,
exp(j2π(n − 1)θk ) for n = 1, . . . , N . We also define
Ψ(θ) , [ψ(θ1 ), · · · , ψ(θdim(θ) )]. The measurement matrix
Φ ∈ CM ×N is either the identity matrix (M = N , complete
data case) or made of a subset of rows of a diagonal matrix
(M < N , incomplete data case). The vector w is a white
Gaussian noise vector with component variance β. The LSE
problem is that of recovering the model order K along with
the frequency θ̃k ∈ [0, 1] and coefficient α̃k ∈ C of each
component k = 1, . . . , K.
A. Estimation Model
The estimation model and inference approach we present
in the following are adaptions of ideas currently available
in the literature. We have carefully combined these ideas to
obtain an iterative scheme which can be implemented with low
complexity as described in Secs. III and IV, while achieving a
performance comparable to that of state-of-the-art algorithms.
Our algorithm is based on Bayesian inference in an estimation model which approximates (1). Specifically, to enable

CN(αk ; 0, γk ),

(5)

k=1

(1)

k=1

K
max
Y

where γk ∈ [0, ∞) is the active-component variance. Sparsitypromoting priors have previously been used for both basis
selection [25] and LSE [15]. The Bernoulli-Gaussian prior
structure that we have adopted above was first introduced in
[26] and used for LSE in [16].
Even though each αk is modelled as Gaussian in (5), the
prior specification is significantly more general than that because the variance of each component is estimated through γk .
In the numerical investigation we demonstrate that our method
works well even when the true density of each coefficient is
not Gaussian.
We finally use an independent and identically distributed
(i.i.d.) uniform prior on the entries in θ:
p(θ) =

K
max
Y
k=1

p(θk ) =

K
max
Y

1 = 1.

(6)

k=1

If further prior information about the frequencies is available,
it can easily be incorporated through p(θ).
1 Since we can never expect to estimate more parameters than the number
of observed observations, we select Kmax = M in our implementation.

HANSEN et al.: SUPERFAST LINE SPECTRAL ESTIMATION

3

B. Approach
By integrating the component coefficients we obtain the
marginal likelihood
Z
p(y|z, θ; β, γ) = p(y|α, z, θ; β)p(α; γ) dα
= CN(y; 0, C)

(7)

with C , βI + A(θz )Γz AH (θz ) and Γz , diag(γz ).
Based on the marginal likelihood we can write the objective
L(z,ζ, β, θ, γ) , − ln p(z, θ|y; β, γ, ζ)

(θẑ ,γẑ )

= ln |C| + y H C−1 y
K
max
X

(zk ln ζ + (1 − zk ) ln(1 − ζ)) + const. (8)

k=1

The variables (z, θ) and model parameters (β, γ, ζ) are estimated by minimizing (8), i.e., we seek the maximum aposteriori (MAP) estimate of (z, θ) and the ML estimate of
(β, γ, ζ). Our algorithm employs a block-coordinate descent
method to find a local minimum (or saddle point) of (8).
For fixed z the first two terms in (8) are equal to the
objective function of stochastic ML [2], and our approach can
therefore be viewed as stochastic ML extended with a variable
model order.
When the above estimates have been computed, the estimated model order is given by the number of active components, i.e. K̂ = ||ẑ||0 , and the entries of θ̂ẑ are the
estimated frequencies. The corresponding coefficients αẑ can
be estimated as follows. First, write the posterior of α as
p(α|y, ẑ, θ̂; β̂, γ̂) ∝ p(y|α, ẑ, θ̂; β̂)p(α; γ̂)
Y
∝ CN(αẑ ; µ̂, Σ̂)
CN(αk ; 0, γ̂k ), (9)
{k:ẑk =0}

where
µ̂ , β̂ −1 Σ̂AH (θ̂ẑ )y

−1
.
Σ̂ , β̂ −1 AH (θ̂ẑ )A(θ̂ẑ ) + Γ̂−1
ẑ

(10)
(11)

As expected the posterior of the coefficients corresponding to
inactive components (those for which ẑk = 0) coincides with
their prior. These are not of interest (they are inconsequential
in the model (2)) and integrating them out gives a Gaussian
posterior over αẑ . If a point estimate of αẑ is needed, the
MAP (which is also the LMMSE) estimate α̂ẑ = µ̂ can be
used2 .

C. Derivation of Update Equations
As mentioned, our algorithm is derived as a blockcoordinate descent method applied on L in (8). The estimates
are updated in the following blocks: ẑ, ζ̂, β̂ and (θ̂ẑ , γ̂ẑ ).
Each update is guaranteed not to increase L. We note that
the frequencies and variances of inactive components (those
2 Note that for computational convenience we write µ̂ = γ̂
ẑ
is defined by (16). See the text following (26).

L(θẑ , γẑ ) = ln |C| + y H C−1 y + const.,
so we need to solve (θ̂ẑ , γ̂ẑ ) = arg min L(θẑ , γẑ ).

= − ln p(y|z, θ; β, γ)p(z; ζ)p(θ) + const.

−

for which ẑk = 0) are not updated, as L does not depend on
these variables.
1) Estimation of frequencies and coefficient variances:
Even when all remaining variables are kept fixed, it is not
tractable to find the global minimizer of L with respect to the
vector of active component frequencies θẑ and variances γẑ .
We therefore resort to a numerical method. Writing only the
terms of (8) which depend on θẑ , we have

q, where q

In [10] a similar optimization problem involving only the
frequencies is solved by Newton’s method. Directly applying
that approach in our case leads to high computational complexity. Methods based on gradient descent have also been
proposed [14], but we have observed that using this approach
leads to slow converge. As we are concerned with computational speed in this paper, we instead use the limited memory
Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm [27].
This algorithm only requires evaluation of the objective function and its gradient. In the following we demonstrate how
these evaluations can be performed with low complexity. At
the same time the per-iteration of L-BFGS is linear in K̂,
namely O(J K̂), where J is the number of saved updates
used in L-BFGS. In our implementation we use J = 10. We
have observed that L-BFGS converges in a small number of
iterations.
The L-BFGS algorithm requires an initial estimate of the
Hessian of L(θẑ , γẑ ), which is subsequently updated in each
iteration of the algorithm. Every update of the activation variable ẑ results in a change in the dimension of the Hessian (the
number of variables in L(θẑ , γẑ ) changes). This means that
the implicit estimate of the Hessian in the L-BFGS algorithm
is reinitialized rather frequently in our estimation scheme. As
a result, the degree of accuracy of the initialization of the
Hessian has a significant impact on the convergence speed of
the algorithm. We therefore propose to initialize L-BFGS with
a diagonal approximation of the Hessian. As shown below,
the diagonal entries of the Hessian can be obtained with low
computational complexity.
The initial estimate of the Hessian must be positive definite.
This is only achieved when all diagonal entries are positive.
Those entries of the diagonal Hessian which are negative
are therefore replaced with the following values: For entries
corresponding to frequency variables we use (50N )2 as the
diagonal Hessian and for the entries corresponding to the
variance of the kth component we use [γ̂ẑ ]−2
k . These heuristic
values have been determined by considering a diagonally
scaled version of the optimization problem (see [28, Sec. 1.3]).
Here follows the required first- and second-order partial
derivatives of L(θẑ , γẑ ) evaluated at the current estimates
(θ̂ẑ , γ̂ẑ ) (see [15] for some hints on how these are obtained):
∂L
= 2[γ̂ẑ ]k Im{tk − qk∗ rk }
∂[θẑ ]k
∂L
= sk − |qk |2
∂[γẑ ]k

(12)
(13)

4

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. XX, NO. XX, MONTH, 2018


∂2L
= 2[γ̂ẑ ]k Re xk − vk + [γ̂ẑ ]k (t2k − xk sk )
2
∂[θẑ ]k
+ [γ̂ẑ ]k (xk |qk |2 + sk |rk |2 − 2tk rk qk∗ )
+ (uk qk∗ − |rk |2 )
∂2L
= 2sk |qk |2 − s2k ,
∂[γẑ ]2k

Z
= − ln

(14)
(15)

q , ΨH (θ̂ẑ )ΦH Ĉ−1 y

(16)

r , ΨH (θ̂ẑ )DΦH Ĉ−1 y


s , diag ΨH (θ̂ẑ )ΦH Ĉ−1 ΦΨ(θ̂ẑ )


t , diag ΨH (θ̂ẑ )DΦH Ĉ−1 ΦΨ(θ̂ẑ )

(17)
(18)
(19)

p(y, αẑ |ẑ, θ̂; β, γ̂)
dαẑ + const.,
f (αẑ )
(24)
R
where f (αẑ ) ≥ 0 is a function which fulfills f (αẑ ) dαẑ =
1. The inequality follows from Jensen’s inequality.
≤−

f (αẑ ) ln

Following EM, we select f (αẑ ) = p(αẑ |y, ẑ, θ̂; β̂ i−1 , γ̂),
where β̂ i−1 denotes the previous noise variance estimate.
Denote the upper bound on the right-hand side of (24) by
Q(β; β̂ i−1 ) and insert f (αẑ ) to get

Q(β; β̂ i−1 ) = M ln β + β −1 tr Σ̂AH (θ̂ẑ )A(θ̂ẑ )
+ β −1 ||y − A(θ̂ẑ )µ̂||2 + const., (25)

(20)
(21)
(22)

The notation diag(·) denotes a vector composed of the diagonal
entries of the (matrix) argument. The matrix Ĉ is that in (7)
evaluated at θ̂ẑ , γ̂ẑ and β̂. We have defined
 the diagonal matrix
D , diag [0, 2π, 4π, . . . , (N − 1)2π]T . In Secs. III and IV
we discuss how the vectors (16)–(22) can be calculated with
low computational complexity.
2) Estimation of activation probability: With all other
variables fixed,
(8) is a convex function of
 2 the objective

∂ L
ζ ∈ [0, 1/2] ∂ζ 2 > 0 . The global minimizer is then found
by differentiating and setting equal to zero. Considering the
constraints on ζ, we update it as


1 ||ẑ||0
.
(23)
ζ̂ = min ,
2 Kmax
3) Estimation of noise variance: Even when keeping all
remaining variables fixed at their current estimate, the globally
minimizing noise variance β in (8) cannot be found in closed
form. An obvious alternative approach would be to incorporate
the estimation of β into L-BFGS together with the estimation
of θẑ and γẑ . However, we have observed this approach to
exhibit slow convergence because the objective function can
be rather “flat” in the variable β (the gradient is small far away
from any stationary point).
In sparse Bayesian learning [25] a similar estimation problem is solved successfully via the expectation-minimization
(EM) algorithm. To use EM, we need to reintroduce α into
the estimation problem. In order to show how EM is integrated
into our coordinate-block descent method and that the update
of β̂ is guaranteed not to increase (8), it is the easiest to
directly use the upper bound associated with EM (see [29]
for a derivation of EM which takes a similar approach).
The updated estimate of β is the minimizer of an upper
bound on the objective function (8). To obtain the upper bound
we write the terms of the objective function which depend on
β, with all other variables kept fixed at their current estimates:
L(β) = − ln p(y|ẑ, θ̂; β, γ̂) + const.

p(y, αẑ |ẑ, θ̂; β, γ̂)
dαẑ + const.
f (αẑ )

Z

where we have defined vectors

u , ΨH (θ̂ẑ )D2 ΦH Ĉ−1 y


v , diag ΨH (θ̂ẑ )D2 ΦH Ĉ−1 ΦΨ(θ̂ẑ )


x , diag ΨH (θ̂ẑ )DΦH Ĉ−1 ΦDΨ(θ̂ẑ ) .

f (αẑ )

where we have used (9) to evaluate expectations involving αẑ
and µ̂ and Σ̂ are calculated from (10)–(11) based on β̂ i−1 . It
is easy to show that the upper bound has a unique minimizer,
which is used as the updated estimate of the noise variance:
!

tr Σ̂AH (θ̂ẑ )A(θ̂ẑ ) + ||y − A(θ̂ẑ )µ̂||2
i
β̂ = max εβ ,
.
M
(26)
To allow low-complexity calculation of β̂ i we use Woodbury’s
matrix inversion identity to show that µ̂ = γ̂ẑ
q and
 PK̂
H
i−1
tr Σ̂A (θ̂ẑ )A(θ̂ẑ ) = k=1 (β̂ sk [γ̂ẑ ]k ).
The update (26) could be applied repeatedly since an
improved upper bound is used each time. Since we have not
observed any advantages by doing so, we simply perform
the update (26) once for each pass in the block-coordinate
descent algorithm. We also note that even though EM is known
to be prone to slow convergence speed, we have observed
empirically that the estimate of β converges fast, typically
within 10 iterations.
It can easily be shown that with the chosen f (αẑ ), the
inequality in (24) holds with equality at β = β̂ i−1 . It then
follows that the new estimate of β does not increase the
value of the objective function (see the proof of Lemma 4
in Appendix A).
4) Deactivation of Components: We now describe the activation and deactivation of components, which is performed
by the single most likely replacement (SMLR) detector [26].
SMLR has previously been demonstrated to perform well for
LSE [14]–[16], [23].
First we write the terms of (8) which depend on the variables
pertaining to the kth component and fix all other variables at
their current estimate. Based on Woodbury’s matrix inversion
identity and the determinant lemma we get (see [24] for
details)
|q∼k (θk )|2
+ s∼k (θk )
!!
1 − ζ̂
+ const., (27)
+ ln (1 + γk s∼k (θk ))
ζ̂

L(zk , θk , γk ) = zk

−

γk−1

HANSEN et al.: SUPERFAST LINE SPECTRAL ESTIMATION

5

with
q∼k (θk ) , ψ H (θk )ΦH Ĉ−1
∼k y
H

s∼k (θk ) , ψ (θk )Φ

H

(28)

Ĉ−1
∼k Φψ(θk ),

where Ĉ∼k , β̂I+A(θ̂ẑ∼k )Γ̂ẑ∼k AH (θ̂ẑ∼k ) and ẑ∼k is equal
to ẑ with the kth entry forced to zero. The matrix Ĉ∼k is thus
the marginal covariance matrix of the observation vector with
the kth component deactivated.
To evaluate if an active component should be deactivated,
we test if the objective L is increased by doing so, i.e., we
test if L(zk = 0, θ̂k , γ̂k ) < L(zk = 1, θ̂k , γ̂k ). This gives the
deactivation criterion for the kth component:
!


|q∼k (θ̂k )|2
1 − ζ̂
− ln 1 + γ̂k s∼k (θ̂k ) < ln
. (29)
γ̂k−1 + s∼k (θ̂k )
ζ̂
This criterion is evaluated for currently active components,
i.e., for k which has corresponding ẑk = 1.
For computational convenience we note that we can obtain
q∼k (θ̂k ) and s∼k (θ̂k ) from q and s with low complexity. First,
write Ĉ∼k = Ĉ − γ̂k Φψ(θ̂k )ψ H (θ̂k )ΦH and use Woodbury’s
identity to obtain
qi
q∼k (θ̂k ) =
1 − γ̂k si
si
s∼k (θ̂k ) =
,
1 − γ̂k si
where qi and si are the ith entries of (16) and (18) with i
denoting the index for which [θ̂ẑ ]i = θ̂k .
5) Component Activation: We now describe a method to
decide if a deactivated component should be activated. This
also involves estimating the frequency and variance of this
component, because no meaningful such estimates are available before the component is activated. Any of the deactivated
components are equally good candidates for activation. In the
following k refers to an arbitrary value for which ẑk = 0. If
no such k exists all components are already activated and the
activation step is not carried out.
Our method is again based on the expression (27). Inspired
by [16], let γ̄ denote the average of the entries in γ̂ẑ . Define
the change in the objective obtained from setting ẑk = 1,
θ̂k = θk , γ̂k = γ̄:
∆L(θk ) = L(1, θk , γ̄) − L(0, θk , γ̄)
= ln (1 + γ̄s∼k (θk ))

1 − ζ̂
ζ̂

(30)
!
−

|q∼k (θk )|2
−1
γ̄ + s∼k (θk )

Note that the last term in (30) does not depend on θk or γ̄.
Then the frequency is found by maximizing the decrease in
the objective, i.e.,
θ̂k = arg min ∆L(θk ),

(31)

θk ∈G

where G is a grid of L equispaced values, i.e., G ,
{0, 1/L, . . . , 1 − 1/L}. The restriction of the estimated frequencies to a grid does not mean that the final frequency
estimates lie on a grid, because they are refined to be in [0, 1)
in subsequent updates of the frequency vector. For this reason,
the choice of L does not have any impact on the estimation

accuracy, provided that it is sufficiently large.3 In Sec. III and
IV we show how q∼k (θk ) and s∼k (θk ) can be evaluated with
low complexity for all θk ∈ G, such that the minimization can
be performed by means of an exhaustive search over G.
The activation procedure continues only if a decrease in the
objective can be obtained by activating a component at θ̂k ,
i.e., if ∆L(θ̂k ) < −εL . The inclusion of the constant εL > 0
is purely technical, as it simplifies our convergence analysis.
It can be chosen arbitrarily small and we select it as machine
precision in our implementation.
After estimating the frequency, the component variance is
selected as γ̂k = arg min ∆L(1, θ̂k , γk ). Using an approach
γk

similar to [24], this minimizer can be shown to be
(
(θk )|2
|q∼k (θk )|2 −s∼k (θk )
> 1,
if |qs∼k
s2∼k (θk )
∼k (θk )
γ̂k =
0
otherwise.

(32)

The component is only activated if4 γ̂k > 0.
It is instructive to explore the activation criterion ∆L(θ̂k ) <
−εL in detail. Since εL is machine precision, we ignore
it (εL = 0) for simplicity. The activation criterion can be
rewritten to the form
!
!

 1 − ζ̂
|q∼k (θ̂k )|2
1
> 1+
ln 1 + γ̄s∼k (θ̂k )
.
s∼k (θ̂k )
γ̄s∼k (θ̂k )
ζ̂
(33)
Denote the left-hand side of (33) as κk . This quantity can be
interpreted as the signal-to-noise ratio of the kth component
[30], [31]. If the sparse Bayesian learning (SBL) model is
used for sparsity promotion, an activation criterion of the from
κk > 1 is obtained [30], [31]. Algorithms using the activation
criterion κk > 1 are known to be prone to the activation
of “artefact” components with very small γ̂k and α̂k at what
seems to be arbitrary frequencies θ̂k . The right-hand side of
(33) is always larger than one and this helps reduce the number
of artefacts which are activated, as demonstrated in [16]. This
favorable phenomenon is caused by the use of the average γ̄
in the definition of ∆L(θk ) (as opposed to inserting γ̂k from
(32), which resembles the SBL approach).
Even still, we have observed the activation of a few artefact
components in our numerical investigations. We therefore
follow the same idea as [30], [31] and heuristically adjust the
criterion (33) to obtain
|q∼k (θk )|2
>
s∼k (θk )


1+

1
γ̄s∼k (θk )


ln (1 + γ̄s∼k (θk ))

1 − ζ̂
ζ̂

!
+ τ,
(34)

where τ ≥ 0 is some adjustment of the threshold. Specifically
we select τ = 5, cf. the numerical study in Sec. VI-B. Our
numerical experiments show that this simple approach is very
effective at avoiding the inclusion of small spurious components. Since the heuristic criterion (34) is stricter than the
3 A numerical investigation (not reported here) shows that the algorithm is
invariant to the choice of L, provided that L ≥ 2N . In our implementation
we use L equal to 8N rounded to the nearest power of 2.
4 When γ̂ = 0 the kth component is effectively deactivated because the
k
corresponding coefficient αk has a zero-mean prior with zero variance, see
(5). The effective deactivation is also seen in the definition of C in (7) and
it further manifests itself as µ̂k = 0 in (10).

6

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. XX, NO. XX, MONTH, 2018

criterion ∆L(θ̂k ) < −εL , it is guaranteed that the activation
of a component decreases the objective function.
D. Outline of the Algorithm and Implementation Details
The algorithm proceeds by repeating the following steps
until convergence:
1) Check if any components can be activated via the
procedure described in Sec. II-C5.
2) Re-estimate the activation probability ζ via (23).
3) Re-estimate the noise variance via (26).
4) Repeat:
4a) Perform a single L-BFGS update of the estimated
vectors of active component frequencies θẑ and
variances γẑ , as described in Sec. II-C1.
4b) Check if any components can be deactivated via
(29).
The algorithm terminates when the change in the objective (8)
between two consecutive iterations is less than M 10−7 .
In step 1) and 4b) the check for component (de)activation is
repeated until no more components can be (de)activated. The
updates in step 4) are iterated until either the approximated
squared Newton decrement of the L-BFGS method is below
M 10−8 or at most 5 times.
The observant reader will have noticed that the minimization
over (θẑ , γẑ ) must be constrained to γk ≥ 0 for all k. It turns
out that this constraint can be handled in a simple manner:
Notice that the deactivation criterion (29) is always fulfilled for
γ̂k sufficiently small. The constraint is therefore never active
at the solution. We therefore simply need to restrict the linesearch performed in L-BFGS such that the no entry in γ̂ẑ ever
becomes negative. If any γ̂k approaches (or becomes equal
to) zero, it is deactivated in step 4b). Note that this approach
resembles that of L-BFGS for box constraints [32], except that
the deactivation of variables for which the constraint is active
happens automatically in our algorithm.
The algorithm is initialized with all components in the
deactivated stage (i.e. ẑ = 0). The initial values of the entries
in θ̂ and γ̂ do not matter, since they are assigned when their
corresponding component is activated (see Sec. II-C5). The
noise variance is initialized to β̂ = 0.01||y||2 /M (1 % of the
energy in y is assumed to be noise). The activation probability
is initialized to ζ̂ = 0.2.
In Appendix A we discuss in detail the convergence properties of our algorithm. The findings are summarized here. We
show that our algorithm terminates in finite time and that the
estimates of z, ζ and β are guaranteed to converge. We denote
the limit points as z̄, ζ̄ and β̄. When these estimates have
converged, our algorithm reduces to a pure L-BFGS scheme
which estimates (θz̄ , γz̄ ). Due to the non-convexity of the
objective function, we cannot guarantee convergence of LBFGS (see [33]). Despite of this, we have never observed nonconvergence of our algorithm. In our experiments it always
converged to a local minimum of the objective function. We
therefore rely on the vast amount of experimental validation
of the convergence of L-BFGS and assume convergence to a
stationary point. In particular, we have the following theorem.

Theorem 1: Assume that L-BFGS in step 4a) converges to
a stationary point of (θz̄ , γz̄ ) 7→ L(z̄, ζ̄, β̄, θz̄ , γz̄ ). Then the
sequence of estimates obtained by our algorithm converges.
Further, the limit point is a stationary point of (ζ, β, θ, γ) 7→
L(z̄, ζ, β, θ, γ), in the sense that the Karush-Kuhn-Tucker
necessary conditions for a minimum are fulfilled.
Proof: See Appendix A.
E. Initial Activation of Components
When the number of sinusoids K in the observed signal (1)
is high, the algorithm spends significant computational effort
activating components (step 1). This is because each time a
component is activated, the values q∼k (θk ) and s∼k (θk ) must
be evaluated for all θk ∈ G to calculate (31). To alleviate
the computational effort of building the initial set of active
components, we propose to let the first few iterations use an
approximate scheme for activating components in place of step
1). The approximate activation scheme proceeds as follows:
1) Calculate q∼k (θ) and s∼k (θ) for all θ ∈ G, where k is
the index of a deactivated component.
2) Evaluate ∆L(θ) (30) for all θ ∈ G.
3) Find the local minimizers of ∆L(θ), i.e., find the values
of θ for which ∆L(θ) ≤ ∆L(θ0 ) with θ0 being any of the
two neighbouring grid-points of θ. The local minimizers
are candidate frequencies.
4) Activate a component at those candidate frequencies for
which the following criteria are fulfilled:
• The component activation criterion (33) is fulfilled.
• The component variance (32) is non-zero.
• The decrease in the objective obeys ∆L(θ) ≤
∆Lmin /5, where ∆Lmin is the largest decrease
obtained from activating a component at another
candidate frequency (in the current iteration).
• All other currently active components have frequency estimates located at least5 0.05N −1 apart
from the candidate frequency.
The above method is a heuristic scheme, which quickly builds
a set of activated components. Typically this set is close to the
final result and only a few (in our setup less than 15 in most
cases) iterations are need before convergence.
III. S UPERFAST M ETHOD (C OMPLETE O BSERVATIONS )
The algorithm presented above has rather large computational complexity, in particular due to the inversion of C
and the calculation of q∼k (θ) s∼k (θ) for all θ ∈ G. In this
section we discuss how all updates of the algorithm can be
evaluated with low computational complexity by exploiting
the inherent structure of the problem. In particular we discuss
how to evaluate ln |C|, y H C−1 y, q, r, s, t, u, v, x and
q∼k (θ), s∼k (θ) for all θ ∈ G.
The method presented here is only applicable when the
complete observation vector is available, i.e., when Φ = I,
M = N and A(θ) = Ψ(θ). In this case the observation
vector y is a wide-sense stationary process and its covariance
5 For the distance measure we use the wrap-around distance on [0, 1) defined
as d(x, y) , min(|x − y|, 1 − |x − y|) for x, y ∈ [0, 1).

HANSEN et al.: SUPERFAST LINE SPECTRAL ESTIMATION

7

matrix C is Hermitian Toeplitz. Low-complexity algorithms
for inverting such matrices are available in the literature. We
also rely on fast Fourier transform (FFT) techniques.
Our approach is based on the Gohberg-Semencul formula
[17], [20], which states that the inverse of the Hermitian
Toeplitz matrix C can be decomposed as

−1
H
H
C−1 = δN
(35)
−1 T1 T1 − T0 T0 ,
where the entries of T0 and T1 are
[T0 ]i,k = ρi−k−1 ,

transform6 (NUFFT) [21], [22]. Then q, r and u are easily
found in time O(N log N ) (assuming the decomposition (35)
has already been calculated).
C. Evaluating s, t, v and x
Turning our attention to s, we follow [19] and note that
(recall that we assume Φ = I)
sk = [ΨH (θ̂ẑ )C−1 Ψ(θ̂ẑ )]k,k
N
−1


X
=
ωs (i) exp j2πi[θ̂ẑ ]k

(38)

i=−(N −1)

[T1 ]i,k = ρN −1+i−k
for i, k = 1, . . . , N . Note that ρi = 0 for i < 0 and
i > N − 1; thus T0 is strictly lower triangular and T1
is unit upper triangular (ρN −1 = 1). The values δi and ρi
for i = 0, . . . , N − 1 can be computed with a generalized
Schur algorithm in time O(N log2 N ) [17]. Alternatively,
the Levinson-Durbin algorithm can also be used to obtain
the decomposition in time O(N 2 ). The latter algorithm is
significantly simpler to implement and is faster for small N .
In [18] it is concluded that the Levinson-Durbin algorithm
requires fewer total operations than the generalized Schur
algorithm for N ≤ 256.

for k = 1, . . . , K̂ where K̂ is the number of entries in θ̂ẑ .
The function ωs (i) gives the sum over the ith diagonal, i.e.,
min(N −1−i,N −1)

X

ωs (i) =

[C−1 ]q+1,q+i+1 .

(39)

q=max(0,−i)

It is obvious that (38) can be calculated for all k = 1, . . . , K̂
via a NUFFT when the values ωs (i) are available.
To evaluate t, v and x we follow a similar approach and
note that the entries of these vectors can be written as (38)
with ωs (i) replaced by
min(N −1−i,N −1)

ωt (i) =
H

A. Evaluating y C

−1

X

[DC−1 ]q+1,q+i+1

(40)

q=max(0,−i)

y and ln |C|

min(N −1−i,N −1)

To calculate the value of the objective function (8) we need
to find y H C−1 y and ln |C|. Inspecting (35) it is clear that
matrix-vector products involving T0 and T1 are convolutions.
These can be implemented using FFT techniques. The product
y H C−1 y can thus be calculated in O(N log N ) time when
{ρi } and δN −1 are known.
The matrix C is Hermitian positive definite and can therefore be factorized uniquely as
C = LBLH ,

N
−1
X

ln δi .

X

[D2 C−1 ]q+1,q+i+1

(41)

[DC−1 D]q+1,q+i+1 ,

(42)

q=max(0,−i)
min(N −1−i,N −1)

ωx (i) =

X
q=max(0,−i)

respectively. In Appendix B we demonstrate how {ωs (i)},
{ωt (i)}, {ωv (i)} and {ωx (i)} can be obtained through length2N FFTs using the decomposition (35).

(36)

with L being unit lower triangular. The diagonal matrix B is
computed with the generalized Schur algorithm. Its diagonal
entries are given by δi for i = 0, . . . , N − 1 [17]. Since the
determinant of a triangular matrix is the product of its diagonal
entries, we have
ln |C| =

ωv (i) =

(37)

i=0

It follows that once the generalized Schur algorithm has been
executed, the objective function (8) can easily be found.

B. Evaluating q, r and u
Note that C−1 y can be evaluated with FFT techniques
using (35). We recognize that matrix-vector products involving
ΨH (θ̂ẑ ) are Fourier transforms evaluated off the equispaced
grid. Such products are approximated to a very high precision in time O(N log N ) using the non-uniform fast Fourier

D. Evaluating q∼k (θ) and s∼k (θ) for all θ ∈ G
To calculate the frequency of the component processed in
the activation stage, q∼k (θ) and s∼k (θ) must be evaluated for
all θ ∈ G, where G is a grid of L equispaced points. Defining
the vector of gridded frequencies θ G , [0, 1/L, . . . , (L −
1)/L]T , we need to find
q G , ΨH (θ G )C−1 y,

sG , diag ΨH (θ G )C−1 Ψ(θ G ) .
We have used the fact that in the beginning of the activation
step the kth component is deactivated and thus C∼k = C.
6 The NUFFT calculates the Fourier transform at arbitrary points (not
lying on an equispaced grid) by interpolation combined with an FFT. It
is an approximation, which can be made arbitrarily accurate by including
more points in the interpolation. The NUFFT achieves a time complexity of
O(N log N + K), where K is the number of off-the-grid frequency points at
which it is evaluated. For K ≤ N this complexity is equal to that of the FFT,
but the constant hidden in the big-O notation is much higher for the NUFFT.
We have found that for N ≥ 512 significant speedups can be achieved by
using the NUFFT over a direct computation of A(θ̂ẑ ) and evaluation of the
matrix-vector products involving this matrix. In particular the speedup arises
from the fact that A(θ̂ẑ ) no longer needs to be formed.

8

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. XX, NO. XX, MONTH, 2018

Since G is an equispaced grid, products with ΨH (θ G ) can
be evaluated as a length-L FFT. The vector q G is therefore
easy to find. Rewriting sG in the form (38), it is seen
that sG can also be evaluated as a length-L FFT. These
computations have time-complexity O(L log L) (assuming the
decomposition (35) has already been calculated).
E. Algorithm Complexity
In summary, the time complexity of each iteration in the
algorithm described in Sec. II is dominated by either the
calculation of {ρi } and δN −1 with the generalized Schur
algorithm or the calculation of q G and sG (we assume K̂ ≤
M = N ≤ L). With our choice L = 8N we have complexity
per iteration of O(N log2 N ).
Also note that all computations involving Ψ(θ̂ẑ ) are performed using the NUFFT. This matrix therefore does not need
to be stored, so our algorithm only uses a modest amount of
memory.
IV. S EMIFAST M ETHOD (I NCOMPLETE O BSERVATIONS )
The method presented in Sec. III is not applicable when an
incomplete observation vector is available, i.e., when Φ 6= I.
In the following we introduce a computational method, which
can be used when Φ is a subsampling and scaling matrix,
i.e., when Φ ∈ CM ×N consists of M rows of a diagonal
matrix.7 With this method we can still obtain an algorithm with
reasonable computational complexity per iteration, assuming
that K̂ is relatively small (a K̂ × K̂ matrix must be inverted).
We coin this algorithm as semifast. For small K̂ the semifast
algorithm is faster than the superfast algorithm of Sec. III and
it may therefore be beneficial to even use it in the complete
data case.
The semifast method is based on the following decomposition of C−1 , obtained using Woodbury’s matrix identity:
C−1 = β̂ −1 I − β̂ −2 A(θ̂ẑ )Σ̂AH (θ̂ẑ )

(43)

with Σ̂ given by (11). We can evaluate Σ̂−1 by noting that
i
h
i
h
AH (θ̂ẑ )A(θ̂ẑ )
= ΨH (θ̂ẑ )ΦH ΦΨ(θ̂ẑ )
i,k

=

M
X

i,k




|Φm,IM (m) |2 exp j2π(IM (m) − 1)(θ̂k − θ̂i ) ,

m=1

(44)
which can be evaluated with a NUFFT in time O(N log N +
K̂ 2 ). Forming Σ̂−1 is then easy and an inversion8 in time
O(K̂ 3 ) is needed to obtain Σ̂. The approach thus hinges on
K̂ being sufficiently small, such that the inverse (really, the
Cholesky decomposition) can be calculated in reasonable time.
7 Let M ⊆ {1, . . . , N } denote the index set of the observed entries and
IM : {1, . . . , M } → M be an indexing. Then Φm,IM (m) , m = 1, . . . , M ,
are the only nonzero elements of Φ.
8 As is customary in numerical linear algebra, we would recommend not
to explicitly evaluate the inverse, but instead use the numerically stabler and
faster approach of calculating the Cholesky decomposition Σ̂−1 = LLH (a
unique Cholesky decomposition exists because Σ̂−1 is Hermitian positive
definite). We need to evaluate matrix-vector products involving Σ̂ which are
easily evaluated from the decomposition by forward-backward substitution.
We can also calculate |Σ̂−1 | directly from the Cholesky decomposition.

A. Evaluating y H C−1 y, ln |C|, q, r and u
Notice that matrix-vector products involving Ψ(θ̂ẑ ) and
ΨH (θ̂ẑ ) can be evaluated using a NUFFT. It then immediately
follows that the values y H C−1 y, q, r and u can be evaluated
using (43) with complexity O(K̂ 2 + N log N ).
To evaluate the objective function (8) we need to calculate
ln |C|. By invoking the matrix determinant lemma we get
X
ln |C| = M ln β̂ +
ln γ̂k + ln |Σ−1 |,
(45)
{k:ẑk =1}

which can be evaluated in time O(K̂) once the Cholesky
decomposition of Σ−1 is known.

B. Evaluating s, t, v and x
As an example, we demonstrate how to evaluate t. We
note that s, v and x can easily be obtained using the same
approach. First, insert (43) into (19) to get


t = β̂ −1 diag ΨH (θ̂ẑ )DΦH A(θ̂ẑ )


− β̂ −2 diag ΨH (θ̂ẑ )DΦH A(θ̂ẑ )Σ̂AH (θ̂ẑ )A(θ̂ẑ ) .
Using the same methodology as for computing Σ̂−1 , the K̂ ×
K̂ matrices AH (θ̂ẑ )A(θ̂ẑ ) and ΨH (θ̂ẑ )DΦH A(θ̂ẑ ) can be
obtained in time O(N log N + K̂ 2 ). Then, t is found by direct
evaluation in time O(K̂ 3 ).

C. Evaluating q∼k (θ) and s∼k (θ) for all θ ∈ G
To calculate the frequency of the component processed in
the activation stage we must evaluate q∼k (θ) and s∼k (θ) for
all θ ∈ G, where G is a grid of L equispaced points. Using
the fact that in the beginning of the activation step the kth
component is deactivated and thus C∼k = C, we obtain the
required quantities by inserting (43) into (16) and (18):


q G = β̂ −1 AH (θ G ) y − A(θ̂ẑ )µ̂

sG = β̂ −1 diag AH (θ G )A(θ G )


− β̂ −2 diag AH (θ G )A(θ̂ẑ )Σ̂AH (θ̂ẑ )A(θ G ) .
It is clear that q G can easily be found using FFT techniques.
note that the first term is a vector with
To obtain sG we first P
M
all entries equal to β̂ −1 m=1 |Φm,IM (m) |2 . The second term
is found by using a NUFFT (see (44)) to form AH (θ̂ẑ )A(θ G ).
Then by using the Cholesky decomposition of Σ̂−1 the second
term can be calculated in time O(LK̂ 2 ).

D. Algorithm Complexity
The above computation is dominated by either the calculation of sG or the length-L FFT involved in calculating q G .
Again with L = 8N we have overall complexity per iteration
O(N K̂ 2 + N log N ).

HANSEN et al.: SUPERFAST LINE SPECTRAL ESTIMATION

9

V. M ULTIPLE M EASUREMENT V ECTORS
The algorithm presented in Sec. II assumes a single measurement vector (SMV). We now discuss an extension to
the case of multiple measurement vectors (MMV) [34]. This
case is of particular importance in array processing where the
number of observation points M is determined by the number
of antennas in the array.9 Typically M is small, which thus
limits estimation accuracy. On the other hand it is often easy
to obtain multiple observation vectors across which the entries
in θ̃ (the true directions of arrivals) are practically unchanged.
The MMV signal model reads
y

(g)

= A(θ̃)α̃

(g)

+w

(g)

,

(46)

where g = 1, . . . , G indexes the observation vectors.
To extend our SMV algorithm to the MMV case we again
impose an estimation model of the form (2) that contains Kmax
components which can be (de)activated based on variables zk ,
k = 1, . . . , Kmax . The likelihood for each of the G observation
vectors then reads
p(y

(g)

|α

(g)

, z, θ; β) = CN(y

(g)

; A(θz )α(g)
z , βI).

(47)

We impose the same prior as used in the SMV case (5) on
each α(g) :
p(α

(g)

; γ) =

K
max
Y

(g)

CN(αk ; 0, γk ).

(48)

k=1

The vectors z and θ are assigned the same priors as in the
SMV case, i.e., as given by (4) and (6). Similarly to the SMV
case, the MMV model has parameters γ, β and ζ.
The objective to be minimized is the marginal likelihood,
which for the MMV model reads
LMMV , − ln

G
Y

p(y (g) |z, θ; β, γ)p(z; ζ)p(θ) + const.,

g=1

=


G 

H
X
ln |C| + y (g) C−1 y (g)
g=1

−

K
max
X

(zk ln ζ + (1 − zk ) ln(1 − ζ)) + const.,

k=1
(g)

where p(y |z, θ; β, γ) = CN(y (g) ; 0, C) with C as in (7).
The posterior probabilities of the coefficient vectors α(g) , g =
1, . . . , G, are given by (9) with y and α replaced by y (g) and
α(g) .
The procedure to estimate the variables θ, z, γ, β and ζ
follows straightforwardly from the method used in the SMV
case. Here we provide a brief discussion of the derivation of
the update equations; refer to Sec. II for details.
To estimate θ and γ the first- and second-order derivatives
of LMMV are needed. Denote the derivative (12) with y
∂L(g)
. Then we have
replaced by y (g) as ∂[θ
ẑ ]k
G

X ∂L(g)
∂LMMV
=
.
∂[θẑ ]k
∂[θẑ ]k
g=1
9 It is worth noting that in array processing the complete data case
corresponds to the very common situation of using a uniform linear array.

A similar result follows for the second-order derivative and
the derivatives with respect to γẑ .
The estimate of ζ is unchanged from the SMV case (23).
To estimate the noise variance β, we write an upper bound
of the same form as (25) and find its minimizer to be


β̂ = max εβ , M −1 tr Σ̂AH (θ̂ẑ )A(θ̂ẑ )
+ (GM )−1

G
X


||y (g) − A(θ̂ẑ )µ̂(g) ||2 ,

g=1

where µ̂(g) is given by (10) with y replaced by y (g) .
To write the activation and deactivation criteria for the
MMV model we rewrite the objective in terms of the parameters of a single component, analogously to (27):
LMMV (zk , θk , γk ) = zk G ln 1 + γk s∼k (θk )
G
X

(g)

|q∼k (θk )|2
1 − ζ̂
+ ln
−
−1
γ
+
s
(θ
)
ζ̂
∼k k
g=1 k



!!
+ const., (49)

(g)

where q∼k (θk ) is given by (28) with y replaced by y (g) . We
omit the details of the activation and deactivation stages, as
they follow straightforwardly from (49) and the description in
Secs. II-C4 and II-C5.
The insightful reader may have noticed that the calculations
required for MMV are very similar to those required for SVM.
In particular, the matrix Ĉ is unchanged and the methods for
calculating matrix-vector products involving Ĉ−1 presented in
Secs. III and IV can be utilized. All expressions involving
y (i.e., q, r, u, q G and y Ĉ−1 y) must be calculated for
each observation vector y (g) . This means that in the case of
complete observations, the generalized Schur algorithm can be
used so that the MMV algorithm has per-iteration complexity
O(N log2 N + GN log N ). With incomplete observations the
semifast method can be used with per-iteration complexity
O(N K̂ 2 + GN log N ).
VI. E XPERIMENTS
A. Setup, Algorithms & Metrics
In our experiments we use the signal model (1). In the
following the wrap-around distance on [0, 1) is used for all
differences of frequencies (see Footnote 5). Unless otherwise
noted, the true frequencies are drawn randomly, such that the
minimum separation between any two frequencies is 2/N .
Specifically, the frequencies are generated sequentially for
k = 1, . . . , K with the kth frequency, θ̃k , drawn from a
uniform distribution on the set {θ ∈ [0, 1) : d(θ, θ̃l ) >
2/N for all l < k}.
The true coefficients in α̃ are generated i.i.d. random, with
each entry drawn as follows. First a circularly-symmetric
complex normal random variable ak with standard deviation
0.8 is drawn. The coefficient is then found as α̃k = ak +
0.2 ej arg(ak ) . The resulting random variable has the property
|α̃k | ≥ 0.2, i.e., all components have significant magnitude.
We use this specification to ensure that all components can
be distinguished from noise. After generating the set of K

frequencies and coefficients, the noise variance β is selected
such that the desired signal-to-noise ratio (SNR) is obtained,
with SNR , ||ΦΨ(θ̃)α̃||2 /(M β).
We compare the superfast LSE algorithm10 with the following reference algorithms: variational Bayesian line spectral
estimation (VALSE) [16]; atomic soft thresholding11 (AST)
[12]; gridless SPICE12 (GLS) [35]; ESPRIT [6], [8]; and a
gridded solution obtained with the least absolute shrinkage
and selection operator (LASSO) solved using SpaRSA13 [36].
The solution to the primal problem of AST [12] directly
provides an estimate of the signal vector h = Ψ(θ̃)α̃. This
solution is known to be biased towards the all-zero solution
(as is also the case with the classical LASSO solution). A
so-called debiased solution can be obtained by recovering the
frequencies from the AST dual and estimating the coefficients
α̃ via least-squares. As in [12], we report here the debiased
solution. If the frequencies are separated by at least 2/N , the
AST algorithm is known to exactly recover the frequencies
in the noise-free case [11]–[13]. In the noisy case no such
recovery guarantee exists, but a bound on the estimation error
of the signal vector h is known [12], [13]. Unfortunately this
error bound does not apply to the debiased solution we report
herein.
We use the variant of GLS [35] which uses SORTE [9]
for model order estimation and MUSIC [7] for frequency
estimation.
ESPRIT requires an estimate of the signal covariance matrix
and of the model order. The former is obtained as the averaged
sample covariance matrix computed from the signal vector
split into N/3 signal vectors of length 2N/N using forwardbackward smoothing. The model order is estimated with
SORTE [9].
The LASSO solution is obtained using a grid of size 8N .
We have observed that no improvement in performance is
achieved with a finer grid. The regularization parameter of
LASSO is selected as proposed in [12] with knowledge of the
true noise variance. We use the debiased solution returned by
the SpaRSA solver.
In the evaluation of the signal reconstruction we have also
included an oracle estimator (denoted Oracle) which obtains
a least squares solution for α̃ with known θ̃.
Three performance metrics are used: normalized meansquared error (NMSE) of the reconstructed signal, block
success rate (BSR) and component success rate (CSR). The
NMSE is defined as
NMSE ,

||Ψ(τ̂ )α̂ − Ψ(τ̃ )α̃||2
.
||Ψ(τ̃ )α̃||2

The BSR is the proportion of Monte Carlo trials in which
the frequency vector θ̃ is successfully recovered. Successful
recovery is understood as correct estimation of the model order
K and that ||d(θ̃, θ̂)||∞ < 0.5/N . The association of the
entries in θ̂ to those in θ̃ is obtained by using the Hungarian
10 We have published our code at github.com/thomaslundgaard/superfast-lse.
It is based on our own implementation of the generalized Schur algorithm and
the NUFFT available at cims.nyu.edu/cmcl/nufft/nufft.html.
11 The code is available online at github.com/badrinarayan/astlinespec.
12 The code has kindly been provided by the authors.
13 The code is available online at lx.it.pt/∼mtf/SpaRSA.

Number of occurences
600
200
400

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. XX, NO. XX, MONTH, 2018

SNR=12dB
SNR=20dB
SNR=30dB

0

10

−5

0

5
10
δ
Fig. 1. Histograms of δ values for different signal-to-noise ratios. Cases where
δ > 0 correspond to cases where an artefact component is activated if the
criterion (33) is used.

method [37] (also known as Munkres assigment algorithm)
minimizing ||d(θ̃, θ̂)||22 .
The BSR can be misleading, since a trial is considered to
be unsuccessful if just a single component is misestimated;
for example if a component is represented in the estimate
by two components with very close frequencies. We therefore
introduce the CSR, defined as follows:
PK̂
PK
k=1 S(θ̂k , θ̃) +
k=1 S(θ̃k , θ̂)
CSR ,
K̂ + K
with the success function S(x, a) , 1[min d(x, ak ) <
k

0.5/N ], where 1[·] denotes the indicator function. The reported
CSR is averaged over a number of Monte Carlo trials. The
CSR takes values in [0, 1]. A CSR of 1 is achieved if, and
only if, all estimated components are in the vicinity of one
or more true components and all true components are in the
vicinity of one or more estimated components.
B. Choosing the Activation Threshold
To determine a sensible value for the activation threshold τ
in (34), the following experiment is conducted. We consider
the complete data case with N = M = 128 and the number of
components is fixed at K = 35, as there is a larger tendency to
activate artefact components for relatively large K/N . The algorithm is provided with the knowledge of Kmax = 35 and the
activation probability is fixed at ζ̂ = 35/128. The algorithm is
run with the activation criterion (33). In this way, the algorithm
in most cases successfully estimates the frequencies without
any artefacts. After the algorithm has terminated, we test if
θ̃ was successfully recovered (as defined above). If so, Kmax
is increased and the procedure for activating a component in
Sec. II-C5 is run and the difference between the left-hand and
right-hand sides of (33) is saved. We refer to this difference
as δ and the criterion (33) can be expressed as δ > 0. In Fig.
1 we show histograms of the value δ obtained from 5, 000
successful recoveries at three different SNR values At each
SNR, the experiment is repeated until the required number
of successful recoveries are obtained; trials without successful
recovery are discarded. Cases where δ > 0 thus correspond
to cases where an artefact would be activated using criterion
(33).
The heuristic criterion (34) corresponds to δ > τ . From
Fig. 1 it is clearly seen that threshold τ = 5 is a sensible
value, which precludes almost all artefact components from
being activated. It is seen that this threshold works well for a

HANSEN et al.: SUPERFAST LINE SPECTRAL ESTIMATION

large range of SNR values. We have not investigated whether
τ = 5 is so large that desired components are precluded
from activation. The results in the following investigations
are all obtained with τ = 5 and the good performance of
our algorithm across this wide selection of scenarios indicates
that the selected τ is not too large.
C. Estimation with Complete Data
In Fig. 2 we show performance results versus SNR. We
first notice that Superfast LSE is on par with or outperforms
all other algorithms in the three metrics shown here for all
SNR values. In the low SNR region no algorithm can reliably
recover the correct model order and the frequencies. In the
plots of the CSR and NMSE, we see that Superfast LSE,
VALSE and AST generally achieve the best approximation of
the true frequencies. There is a small performance gap in terms
of NMSE between Oracle and all other evaluated algorithms
due to the uncertainty in frequency estimation (Oracle knows
the true frequencies).
ESPRIT and GLS are observed to have the weakest performance at low SNR, especially in terms of CSR and NMSE.
Both algorithms use SORTE to estimate the model order
from the eigenvalues of the signal covariance. At low SNR
it is hard to distinguish the signal eigenvalues from the
noise eigenvalues, leading to the observed deterioration in
performance.
At medium to high SNR, BSR of AST is about 0.75. The
algorithm tends to slightly overestimate the model order (not
shown here). We hypothesise that such systematic overestimation of the model order can be avoided by adjusting the regularization parameter used in AST. Doing so would, however,
mean that AST would perform worse in other scenarios. This
is exactly the weakness of methods involving regularization
parameters.
Finally note that LASSO is never able to successfully
estimate the model order, due to the use of a grid. In particular
each true frequency component is estimated by a few nonzero entries at neighbouring gridpoints. It is visible in the
CSR that LASSO indeed estimates frequencies which lie in
the vicinity of the true frequencies. In some applications,
e.g. channel estimation in wireless communications, it is the
reconstructed signal and not the frequencies themself which
are of interest. In this case LASSO may be preferable because
of its simplicity. Due to the grid approximation, LASSO
performs a little worse than the best gridless algorithms in
terms of NMSE.
D. Super Resolution
The ability to separate components beyond the Rayleigh
limit of 1/N is known as super resolution. The results in
Fig. 3 illustrate the super resolution ability of the algorithms.
In this experiment we generate 5 pairs of frequencies with
varying distance between the two frequencies in each pair.
The pairs are well separated (at least 2/N separation between
frequencies which are not in the same pair).
The NMSE performance of Superfast LSE, VALSE and
LASSO is only slightly worse at low separation when compared to the performance at large separation. It is evident that

11

the model order and the frequencies cannot be recovered in
every case (BSR below 1) when the separation is less than
1/N . Since the CSR is close to 1 and the NMSE is close to
that of Oracle, these three algorithms handle closely located
components well, in the sense that a good approximation of
the frequencies is obtained.
AST, GLS and ESPRIT give a CSR below 1 and a rather
large NMSE when the separation is small. This is despite the
fact that GLS and ESPRIT do not show a significantly worse
BSR compared to Superfast LSE. We have observed that this is
because these algorithms significantly underestimate the model
order in some cases, resulting in a large contribution to NMSE.
ESPRIT shows worse super resolution ability than Superfast
LSE, VALSE and GLS (lower BSR for separation below
0.7/N ). This is because a covariance matrix of size 2N/3
is formed, thus reducing the effective signal length.
E. Estimation with Incomplete Data
Fig. 4 reports the performance in the incomplete data
case. The measurement matrix Φ is generated by randomly
selecting M rows of the N × N identity matrix. The set
of observation indices is chosen to include the first and last
indices, while the remaining M − 2 indices are obtained by
uniform random sampling without replacement. Only a subset
of the algorithms are applicable in this case. Our proposed
algorithm is implemented using the techniques described in
Sec. IV. We refer to it as Semifast LSE.
Semifast LSE and VALSE largely show the same performance, while GLS has a slightly higher NMSE for M/N ≤
0.5. The higher NMSE is caused by a few outliers (less
than 1% of the Monte Carlo trials) where GLS significantly
underestimates the model order. LASSO is again observed
to have reasonable NMSE and CSR while being unable to
correctly estimate the set of frequencies (i.e., BSR= 0).
F. Phase Transitions
Inspired by the concept of phase transitions in compressed
sensing, we perform an experiment which shows a similar phenomenon for LSE. In particular we demonstrate that for each
algorithm there is a region in the space of system parameters
where it can almost perfectly recover the frequencies and a
region where it cannot, with a fairly sharp transition between
the two. The results, in terms of BSR, are seen in Fig. 5.
We first note that AST has rather poor performance, which
is consistent with the observation in Fig. 2 that its BSR is
significantly below 1. Turning our attention to VALSE, GLS
and ESPRIT, we see that these algorithms generally do not
deal well with a large number of components, in the sense
that the BSR is significantly below 1 for K/N ≥ 0.15. It is
seen in Fig. 5f that Superfast LSE has the largest region with
high probability of successful recovery (BSR ≥ 0.75).
G. Computation Times
In Fig. 6 and 7 we show algorithm runtimes for varying
problem sizes. Our proposed method uses the superfast and
semifast implementations in Secs. III and IV, respectively. The

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. XX, NO. XX, MONTH, 2018

0
NMSE (dB)
−20

9 12 15 20 25 30
−5
0 3 6 9 12 15 20 25 30
−5
0 3 6 9 12 15 20 25 30
SNR (dB)
SNR (dB)
SNR (dB)
Fig. 2. Simulation results for varying SNR with complete data (Φ = I). The signal length is N = M = 128 and the number of components is K = 10.
Results are averaged over 500 Monte Carlo trials. The legend applies to all plots. Only the NMSE of Oracle is shown.

75
1.

50
1.

25
1.

0.
20
0.
40
0.
55
0.
70
0.
85
1.
00

0

1.
75

1.
5

1.
25

75
1.

50
1.

25
1.

0. 0
20
0.
40
0.
55
0.
70
0.
85
1.
00

Superfast LSE
VALSE
AST
GLS
ESPRIT
LASSO
Oracle

NMSE (dB)
−30 −20 −10

6

Component success rate
0.8
0.9
1

3

0.
20
0.
40
0.
55
0.
70
0.
85
1.
00

0

Block success rate
0.2 0.4 0.6 0.8 1

−5

Superfast LSE
VALSE
AST
GLS
ESPRIT
LASSO
Oracle

−40

0

Component success rate
0 0.2 0.4 0.6 0.8 1

Block success rate
0.2 0.4 0.6 0.8 1

12

0
Semifast LSE
VALSE
GLS
LASSO
Oracle

NMSE (dB)
−30 −20 −10

Component success rate
0 0.2 0.4 0.6 0.8 1

0

Block success rate
0.2 0.4 0.6 0.8 1

Separation (scaled by 1/N )
Separation (scaled by 1/N )
Separation (scaled by 1/N )
Fig. 3. Simulation results for closely23 located components with complete data (Φ = I). The frequencies are generated as 5 pairs (i.e. K = 10) such that
each pair has varying (deterministic) intra-pair separation, while the location of the pairs are generated randomly with non-paired frequencies at least 2/N
apart (i.e., the location of the pairs are generated using a procedure similar to the one described in Sec. VI-A). The signal length is N = M = 128 and the
SNR is 20 dB. Results are averaged over 500 Monte Carlo trials. The legend applies to all plots. Only the NMSE of Oracle is shown.

0.2 0.3 0.4 0.5 0.6
0.8
1
0.2 0.3 0.4 0.5 0.6
0.8
1
0.2 0.3 0.4 0.5 0.6
0.8
1
M/N
M/N
M/N
Fig. 4. Simulation results with incomplete data, i.e., Φ contains M rows of I selected at random. The signal length is N = 128, the SNR is 20 dB and the
number of components is K = 10. Results are averaged over 500 Monte Carlo trials. The legend applies to all plots. Only the NMSE of Oracle is shown.

results are obtained using MATLAB on a 2011 MacBook Pro.
To avoid differences in results originating from the amount
of parallelism achieved in each implementation, MATLAB is
restricted to only use a single computational thread. The part
of the code where each algorithm spends significant time is
implemented as native code via MATLAB’s codegen feature.
For varying N (Fig. 6), we first observe that for small to
moderate problem sizes (N ≤ 210 ) the difference between
LASSO and our proposed algorithms is small (less than 1
second). This difference is mainly due to implementation
details. In the large-N region, Superfast and Semifast LSE
are approximately an order of magnitude faster than LASSO.
We observe that the asymptotic per-iteration complexity of
Superfast and Semifast LSE describes the scaling of the total
runtime well for N ≥ 212 , because the number of iterations
(not shown) stays practically constant. The state-of-the-art
LSE methods VALSE, AST and GLS all have O(N 3 ) time
complexity or worse, which results in very large runtimes even
when the problem size is moderate (e.g. > 100 s for AST and
GLS at N = 512). For large problem sizes, the O(N 3 ) time
complexity of ESPRIT is evident and Superfast/Semifast LSE
and LASSO significantly outperform ESPRIT.

In Fig. 7 we show results illustrating how the computation
time scales with K. In this analysis we assume K̂ = O(K).
First we note that the runtime of LASSO is practically constant
with K. In the complete data case, the per-iteration complexity
of Superfast LSE scales linearly with K. In practice we see a
slower scaling with K, which means that the values of K we
use here are not large enough to reach the asymptotic region.
Simulations with large K cannot be run, because we need
approximately K < N/4 for K̂ = O(K) to hold (cf. Fig. 5).
In the incomplete data case, we see that the runtime of
Semifast LSE increases quickly with K, such that for K >
128 LASSO is faster than Semifast LSE. We do, however,
see that the asymptotic complexity of O(K 3 ) is not reached
in our experiment, because the runtime is dominated by the
calculation of sG , which has complexity O(LK 2 ).
VII. C ONCLUSIONS
We have presented a low-complexity algorithm for line
spectral estimation. Computational methods for both the complete and incomplete data cases have been presented, along
with an extension to the case of multiple measurement vectors.

HANSEN et al.: SUPERFAST LINE SPECTRAL ESTIMATION

13

(b) VALSE

0.6
0.4

Separation (scaled by 1/N )

0.0

1.7

1.5

1.3

1.1

0.7

0.9

1.7

1.5

1.3

1.1

0.9

0.7

0.5

0.3

(f) Contour at 75 % success rate
0.39
Superfast LSE
0.33
VALSE
AST
0.27
GLS
0.21
ESPRIT
0.15
0.09
0.03
0.1

K/N
1.7

1.5

1.3

1.1

0.7

0.9

0.5

0.1

0.2
0.3

K/N
1.7

1.5

1.3

1.1

0.9

0.7

0.5
0.5

0.1

Separation (scaled by 1/N )

0.3

K/N

0.39
0.33
0.27
0.21
0.15
0.09
0.03

1.7

1.5

1.3

1.1

0.9

0.7

0.5

0.3

0.8

Separation (scaled by 1/N )

(e) ESPRIT

0.39
0.33
0.27
0.21
0.15
0.09
0.03

1.0

0.39
0.33
0.27
0.21
0.15
0.09
0.03

Separation (scaled by 1/N )

(d) GLS

0.1

0.3

0.1

K/N
Separation (scaled by 1/N )

K/N

(c) AST

0.39
0.33
0.27
0.21
0.15
0.09
0.03

1.7

1.3

1.5

1.1

0.9

0.7

0.5

0.3

0.1

K/N

(a) Superfast LSE
0.39
0.33
0.27
0.21
0.15
0.09
0.03

Separation (scaled by 1/N )

Fig. 5. Simulation results showing phase transitions with complete data (Φ = I). The plots show block success rate. The set of frequencies are generated as
closely located pairs, following the methodology described in the caption of Fig. 3. The number of pairs are selected such that the specified ratio K/N is
achieved as closely as possible. The signal length is N = M = 128 and the SNR is 20 dB. Results are averaged over 120 Monte Carlo trials.
103

Superfast LSE
LASSO
ESPRIT
VALSE
AST
GLS
O(N log2 N )

102
101
100
10−1
26

212

214

216

218

101
100
10−1
26

220

103

Superfast LSE
LASSO
O(K)

102

Runtime in seconds

Runtime in seconds

210

Semifast LSE
LASSO
VALSE
GLS
O(N log N )

102

28
210
212
214
216
218
220
N
N
Fig. 6. Runtimes in seconds versus problem size N . We show for both complete (left) and incomplete (right) data case. The number of components is
K = 15 and the SNR is 20 dB. Values are averaged over 20 Monte Carlo trials. In the incomplete data case we set M = 0.75N . We also plot (solid black)
the asymptotic per-iteration complexity of Superfast and Semifast LSE.
103

28

Runtime in seconds

Runtime in seconds

103

101
100

102

Semifast LSE
LASSO
O(K 3 )

101
100

10−1
4
8
16
32
64
128
256
512 1,024
64
128
256
512 1,024
K
K
Fig. 7. Runtimes in seconds versus number of components K. We show for both complete (left) and incomplete (right) data case. The problem size is
N = 4096 and the SNR is 20 dB. Values are averaged over 20 Monte Carlo trials. We only show results for algorithms which has runtime lower than 10
s at K = 4. In the incomplete data case we set M = 0.75N = 3072. We also plot (solid black) the asymptotic per-iteration complexity of Superfast and
Semifast LSE.
10−1

4

8

16

32

The proposed algorithm falls in the category of Bayesian
methods for line spectral estimation. Bayesian methods are
widely accepted due to their high estimation accuracy, but
a drawback of this class of methods has historically been
their large computational complexity. In that respect, this work
makes an important contribution in making Bayesian methods
more viable in practice.
At the core of the computational method for the complete
data case lies the application of the Gohberg-Semencul formula to the Toeplitz signal covariance matrix. Many methods
for line spectral estimation have Toeplitz covariance matrices
at their core and we conjecture that the computational com-

plexity of some of them can be drastically reduced by applying
the techniques we have demonstrated in this paper.

Our numerical experiments show that our Superfast LSE
algorithm has very high estimation accuracy. For example,
in Fig. 5 we see that Superfast LSE attains high frequency
recovery rates for a much larger set of scenarios than any
of the reference algorithms. At the same time our algorithm
has so low computation time that it makes highly-accurate
LSE feasible for problems with size much larger than methods
currently available in the literature can practically deal with.

14

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. XX, NO. XX, MONTH, 2018

A PPENDIX A
C ONVERGENCE A NALYSIS
We now discuss the convergence of our proposed blockcoordinate descent algorithm. To do so we introduce an
iteration index i on all estimated variables. Our algorithm
then produces sequences of blocks of estimates denoted as
{ẑ i }, {ζ̂ i }, {β̂ i } and {(θ̂ i , γ̂ i )}. We denote the value of
the objective function at the end of the ith iteration as
Li , L(ẑ i , ζ̂ i , β̂ i , θ̂ i , γ̂ i ). We first note that all updates of
the estimates are guaranteed not to increase the objective, so
{Li } is a non-increasing sequence. Since β ≥ εβ it is also
bounded below and thus it converges. Therefore, our proposed
algorithm terminates in finite time.
Unfortunately, convergence of {Li } does not imply convergence of the sequences of estimates. Even with exact
minimization in each block of variables, block-coordinate
descent on non-convex functions can get stuck in an infinite
cycle [38]. This is further complicated by the fact that our
algorithm only approximately solves the minimization in some
of the blocks.
Proposition 5 in [39] shows that with exact minimization in
each block, block-coordinate descent converges if the objective
function is strictly quasiconvex in all but 2 blocks. The
objective function L is strictly quasiconvex in ζ and β. There is
thus hope that we can prove convergence of our scheme which,
in lieu of computing the exact minimizer, merely has a descent
property in each block. Our approach to show convergence
shares the same overall idea as that in [39], while many details
differ.
To discuss the convergence properties, we first derive a
number of lemmas. Theorem 1 is proved at the end of the
appendix. For notational simplicity we take the convention
that for each i the block-coordinate descent algorithm cycles
through the block updates in the following order: ẑ i , ζ̂ i , β̂ i and
finally (θ̂ i , γ̂ i ), such that for example β̂ i is found based on
(ẑ i , ζ̂ i , β̂ i−1 , θ̂ i−1 , γ̂ i−1 ). This is strictly speaking not how
we have defined our algorithm, but that does not affect the
correctness of our analysis.
Lemma 1: The sequence of estimates has at least one
convergent subsequence, i.e., at least one limit point.
Proof: All variables but γ and β are defined to be in
a closed and bounded set. Since limβ→∞ L = ∞ we can
restrict β to a closed and bounded set determined by the (finite)
initial value of the objective function. A similar argument
holds for each γk . The lemma then follows from the BolzanoWeierstrass theorem.
Lemma 2: The sequence {ẑ i } converges.
Proof: Each activation of a component gives a decrease
in L of at least εL . Since {Li } is lower bounded, there can
only be finitely many activations. Since there cannot be more
deactivations than activations, also the number of deactivations
is finite. There is thus a finite number of changes to ẑ and {ẑ i }
converges. We denote the limit point as z̄.
Lemma 3: The sequence {ζ̂ i } converges. Further, the limit
point ζ̄ is the unique global minimizer of ζ 7→ L(z̄, ζ, β, θ, γ)
for any β, θ and γ.
Proof: The first statement follows from Lemma 2 since
ζ̂ i (23) is only a function of ẑ i . The second statement results

from the fact that ζ̂ i is defined as the global minimizer of
ζ 7→ L(ẑ i , ζ, β̂ i−1 , θ̂ i−1 , γ̂ i−1 ) and this global minimizer does
not depend on β̂ i−1 , θ̂ i−1 and γ̂ i−1 .
Lemma 4: The sequence {β̂ i } converges to the limit point
β̄. Further, for every limit point (z̄, ζ̄, θ̄, γ̄) of the remaining
variables, the limit point β̄ is a local minimum at the boundary
εβ or a stationary point of β 7→ L(z̄, ζ̄, β, θ̄, γ̄).
Proof: To perform this proof we expand our previous notation and denote the upper bound (25) as
Q(β; ẑ i , β̂ i−1 , θ̂ i−1 , γ̂ i−1 ). Then,
L(ẑ i , ζ̂ i , β̂ i−1 , θ̂ i−1 , γ̂ i−1 ) = Q(β̂ i−1 ; ẑ i , β̂ i−1 , θ̂ i−1 , γ̂ i−1 )
≥ Q(β̂ i ; ẑ i , β̂ i−1 , θ̂ i−1 , γ̂ i−1 ) ≥ L(ẑ i , ζ̂ i , β̂ i , θ̂ i−1 , γ̂ i−1 ).
Recalling that {Li } converges, we have
lim L(ẑ i , ζ̂ i , β̂ i−1 , θ̂ i−1 , γ̂ i−1 )

i→∞

− L(ẑ i , ζ̂ i , β̂ i , θ̂ i−1 , γ̂ i−1 ) = 0,
and thus
lim Q(β̂ i−1 ; ẑ i , β̂ i−1 , θ̂ i−1 , γ̂ i−1 )

i→∞

− Q(β̂ i ; ẑ i , β̂ i−1 , θ̂ i−1 , γ̂ i−1 ) = 0. (50)
Reasoning by contradiction, assume that the sequence of
estimates {β̂ i } has two limit points β̄1 and β̄2 , such that
β̄1 6= β̄2 . Let (θ̄, γ̄) be any limit point of {θ̂ i , γ̂ i } (such a
limit point exists due to Lemma 1). Then by (50) we must
have
Q(β̄1 ; z̄, β̄1 , θ̄, γ̄) = Q(β̄2 ; z̄, β̄1 , θ̄, γ̄).

(51)

Recalling the definition of β̂ i , we have that β̄2 uniquely
minimizes Q. Then, since we assumed β̄1 6= β̄2 , we get
Q(β̄1 ; z̄, β̄1 , θ̄, γ̄) > Q(β̄2 ; z̄, β̄1 , θ̄, γ̄),
which contradicts (51). So {β̂ i } has only a single limit point
which we denote as β̄.
To prove the second statement, use (24) to show that
∂
Q(β; z̄, β̄, θ̄, γ̄)
∂β

=
β=β̄

∂
L(z̄, ζ̄, β, θ̄, γ̄)
∂β

.
β=β̄

If β̄ = εβ , we have that the derivatives of Q and thus of L
are positive. It follows that β̄ is a local or global minimum.
If β̄ 6= εβ it is, by definition, a stationary point of Q. It is
therefore also a stationary point of β 7→ L(z̄, ζ̄, β, θ̄, γ̄).
We can now give a proof of the main convergence result.
Proof of Theorem 1: Convergence to a unique limit
follows immediately from the assumption and Lemmas 2, 3
and 4.
To prove the second statement, we first note that L is
constant with respect to those entries of θ and γ for which
∂L
z̄k = 0. It then follows from the assumption that ∂θ
= 0 and
k
∂L
=
0
for
all
k
=
1,
.
.
.
,
K
at
the
limit
point.
Similarly
max
∂γk
from Lemma 3 we have ∂L
∂ζ = 0 at the limit point. If β̄ 6= εβ
we have ∂L
∂β = 0 at the limit point and the result follows
immediately.
If β̄ = εβ the result can be obtained by introducing a

HANSEN et al.: SUPERFAST LINE SPECTRAL ESTIMATION

15

Lagrange multiplier such that the limit point satisfies the
Karush-Kuhn-Tucker conditions.

A PPENDIX B
E FFICIENT E VALUATION OF ωs (i), ωt (i), ωv (i) AND ωx (i)
We derive a low-complexity computation of ωs (i) by first
inserting (35) into (39) to get
ωs (i) =

−1
δN
−1

min(N −1,N −1−i) q
X
X

−ρq−r−1 ρ∗q+i−r−1

r=0

q=max(0,−i)

!
+ ρ∗N −1+r−q ρN −1+r−q−i

(52)

for i = −(N − 1), . . . , N − 1. Then note that since C is
Hermitian we have ωs (−i) = ωs∗ (i). We therefore restrict our
attention to i ≥ 0 in the following. We need the identity
q
N
−1 X
X

zq,r =

−1−q
N
−1 NX
X
q=0

q=0 r=0

zq+k,k ,

(53)

k=0

from which we get (recall that ρi = 0 for i < 0 and i > N −1)
−1
ωs (i) = δN
−1

NX
−1−i

(N − i − q)(ρ∗N −1−q ρN −1−q−i − ρq−1 ρ∗q−1+i ).

q=0

Substituting q = N − 1 − q̄ − i in the first term and q = q̄ + 1
in the second term we finally obtain
ωs (i) =

−1
δN
−1

N
−1
X

cs (q̄, i)ρq̄ ρ∗q̄+i ,

(54)

q̄=0

where cs (q̄, i) , (2 − N + i + 2q̄). The above expression
can be calculated as the sum of two cross-correlations in time
O(N log N ) by using FFT techniques.
To evaluate ωt (i) (40) we again insert (35) and get
min(N −1,N −1−i)

X

−1
ωt (i) = δN
−1

2πq

q=max(0,−i)

q
X

−ρq−r−1 ρ∗q+i−r−1

r=0

!
+

ρ∗N −1+r−q ρN −1+r−q−i

(55)

for i = −(N − 1), . . . , N − 1. Be aware that we do not
have ωt (i) = ωt∗ (−i). Applying (53), performing the same
substitutions as above and following tedious, but straightforward, algebra we finally get
ωt (i) =

2π

N
−1
X

δN −1

q̄=0

ct (q̄, i)ρq̄ ρ∗q̄+i

(56)

with



3+i
ct (q̄, i) , −q̄(q̄ + i) + i N −
2


N −2
+ q̄ 2 + (N − 1) q̄ −
, (57)
2
which again can be evaluated using FFT techniques.

Omitting details, we use a similar approach to find
ωv (i) =

N −1
4π 2 X
cv (q̄, i)ρq̄ ρ∗q̄+i
δN −1 q̄=0

(58)

ωx (i) =

N −1
4π 2 X
cx (q̄, i)ρq̄ ρ∗q̄+i ,
δN −1 q̄=0

(59)

where ωv (−i) = ωv∗ (i). The expression giving ωv (i) is valid
for i ≥ 0, while that giving ωx (i) is valid for i = −(N −
1), . . . , N − 1. We have also defined
cv (q̄, i) = q̄(q̄ + i)2 + (3q̄ − 2N q̄ − q̄ 2 )(q̄ + i)


2 3
7
2
2
+ q̄ + (N − 1)q̄ + N − 3N +
q̄
3
3


3
1
13
− M i (i − N ) + 1
+ (i − N )2 + (i3 − N 3 ) +
2
3
6
cx (q̄, i) = (q̄ 2 + 2q̄ − N q̄)(q̄ + i)


1 3
7
1
2
− q̄ + N − 3N +
q̄ − i3
3
3
6


2
3N − 9N + 7
3
13
1
+
i − N 3 + N 2 − N + 1.
6
3
2
6
R EFERENCES
[1] D. Malioutov, M. Cetin, and A. S. Willsky, “A sparse signal reconstruction perspective for source localization with sensor arrays,” IEEE Trans.
Signal Process., vol. 53, pp. 3010–3022, Aug 2005.
[2] B. Ottersten, M. Viberg, and T. Kailath, “Analysis of subspace fitting
and ML techniques for parameter estimation from sensor array data,”
IEEE Trans. Signal Process., vol. 40, pp. 590–600, Mar. 1992.
[3] R. Carriere and R. L. Moses, “High resolution radar target modeling
using a modified Prony estimator,” IEEE Trans. Antennas Propag.,
vol. 40, pp. 13–18, Jan. 1992.
[4] W. Bajwa, A. Sayeed, and R. Nowak, “Compressed channel sensing:
A new approach to estimating sparse multipath channels,” Proc. IEEE,
vol. 98, pp. 1058–1076, Jun. 2010.
[5] X. Andrade, J. N. Sanders, and A. Aspuru-Guzik, “Application of
compressed sensing to the simulation of atomic systems,” Proc. Nat.
Academy of Sciences, vol. 109, pp. 13 928–13 933, Jul. 2012.
[6] S.-Y. Kung, K. S. Arun, and B. D. Rao, “State-space and singular-value
decomposition-based approximation methods for the harmonic retrieval
problem,” J. of the Optical Soc. of America, vol. 73, pp. 1799–1811,
Dec. 1983.
[7] R. Schmidt, “Multiple emitter location and signal parameter estimation,”
IEEE Trans. Antennas Propag., vol. 34, pp. 276–280, Mar. 1986.
[8] R. Roy and T. Kailath, “ESPRIT - estimation of signal parameters via
rotational invariance techniques,” IEEE Trans. Acoust., Speech, Signal
Process., vol. 37, pp. 984–995, Jul. 1989.
[9] Z. He, A. Cichocki, S. Xie, and K. Choi, “Detecting the number of
clusters in n-way probabilistic clustering,” IEEE Trans. Pattern Anal.
Mach. Intell., vol. 32, pp. 2006–2021, Nov 2010.
[10] L. Hu, Z. Shi, J. Zhou, and Q. Fu, “Compressed sensing of complex
sinusoids: An approach based on dictionary refinement,” IEEE Trans.
Signal Process., vol. 60, pp. 3809–3822, Jul. 2012.
[11] G. Tang, B. N. Bhaskar, P. Shah, and B. Recht, “Compressed sensing
off the grid,” IEEE Trans. Inf. Theory, vol. 59, pp. 7465–7490, Nov.
2013.
[12] B. N. Bhaskar, G. Tang, and B. Recht, “Atomic norm denoising with
applications to line spectral estimation,” IEEE Trans. Signal Process.,
vol. 61, pp. 5987–5999, Dec. 2013.
[13] E. J. Candès and C. Fernandez-Granda, “Super-resolution from noisy
data,” J. of Fourier Anal. and Applicat., vol. 19, pp. 1229–1254, Dec.
2013.
[14] F. Dublanchet, J. Idier, and P. Duwaut, “Direction-of-arrival and frequency estimation using Poisson-Gaussian modeling,” in Proc. IEEE
Int. Conf. Acoust., Speech, and Signal Process., vol. 5, Apr. 1997, pp.
3501–3504.

16

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. XX, NO. XX, MONTH, 2018

[15] T. L. Hansen, M.-A. Badiu, B. H. Fleury, and B. D. Rao, “A sparse
Bayesian learning algorithm with dictionary parameter estimation,”
in Proc. 8th IEEE Sensor Array and Multichannel Signal Process.
Workshop, Jun. 2014, pp. 385–388.
[16] M. A. Badiu, T. L. Hansen, and B. H. Fleury, “Variational Bayesian
inference of line spectra,” IEEE Trans. Signal Process., vol. 65, pp.
2247–2261, May 2017.
[17] G. S. Ammar and W. B. Gragg, “The generalized Schur algorithm for
the superfast solution of Toeplitz systems,” in Rational approximation
and its applications in mathematics and physics, 1987, pp. 315–330.
[18] ——, “Numerical experience with a superfast real Toeplitz solver,”
Linear Algebra and its Applicat., vol. 121, pp. 185–206, Aug. 1989.
[19] B. R. Musicus, “Fast MLM power spectrum estimation from uniformly
spaced correlations,” IEEE Trans. Acoust., Speech, Signal Process.,
vol. 33, pp. 1333–1335, Oct. 1985.
[20] I. Gohberg and I. A. Feldman, Convolution Equations and Projection
Methods for Their Solution, ser. Translations of mathematical monographs. American Mat. Soc., 2005, vol. 41.
[21] L. Greengard and J.-Y. Lee, “Accelerating the nonuniform fast Fourier
transform,” SIAM Review, p. 443–454, 2004.
[22] J.-Y. Lee and L. Greengard, “The nonuniform FFT of Type 3 and its
applications,” J. Comput. Phys., pp. 1–5, 2005.
[23] F. Champagnat, Y. Goussard, and J. Idier, “Unsupervised deconvolution
of sparse spike trains using stochastic approximation,” IEEE Trans.
Signal Process., vol. 44, pp. 2988–2998, Dec. 1996.
[24] M. E. Tipping and A. Faul, “Fast marginal likelihood maximisation for
sparse Bayesian models,” in Proc. 9th Int. Workshop Artificial Intell.
and Stat., Jan. 2003, pp. 3–6.
[25] D. P. Wipf and B. D. Rao, “Sparse Bayesian learning for basis selection,”
IEEE Trans. Signal Process., vol. 52, pp. 2153–2164, Aug. 2004.
[26] J. J. Kormylo and J. Mendel, “Maximum likelihood detection and
estimation of Bernoulli-Gaussian processes,” IEEE Trans. Inf. Theory,
vol. 28, pp. 482–488, May 1982.
[27] J. Nocedal and S. J. Wright, Numerical optimization, 2nd ed. Springer,
2006.
[28] D. P. Bertsekas, Nonlinear programming. Athena scientific Belmont,
1999.
[29] F. Dellaert, “The expectation maximization algorithm,” Georgia Institute
of Technology, Tech. Rep. GIT-GVU-02-20, Feb. 2002.
[30] D. Shutin, B. H. Fleury, and N. Schneckenburger, “Artifact suppression
for super-resolution sparse Bayesian learning,” 2018, submitted to IEEE
Trans. Signal Process.
[31] D. Shutin, T. Buchgraber, S. R. Kulkarni, and H. V. Poor, “Fast variational sparse Bayesian learning with automatic relevance determination
for superimposed signals,” IEEE Trans. Signal Process., vol. 59, pp.
6257–6261, Dec. 2011.
[32] R. H. Byrd, P. Lu, J. Nocedal, and C. Zhu, “A limited memory algorithm for bound constrained optimization,” SIAM Journal on Scientific
Computing, vol. 16, pp. 1190–1208, 1995.
[33] W. F. Mascarenhas, “The BFGS method with exact line searches fails
for non-convex objective functions,” Math. Programming, vol. 99, pp.
49–61, Jan. 2004.
[34] S. F. Cotter, B. D. Rao, K. Engan, and K. Kreutz-Delgado, “Sparse
solutions to linear inverse problems with multiple measurement vectors,”
IEEE Trans. Signal Process., vol. 53, pp. 2477–2488, Jul. 2005.
[35] Z. Yang and L. Xie, “On gridless sparse methods for line spectral
estimation from complete and incomplete data,” IEEE Trans. Signal
Process., vol. 63, pp. 3139–3153, Jun. 2015.
[36] S. J. Wright, R. D. Nowak, and M. A. Figueiredo, “Sparse reconstruction
by separable approximation,” IEEE Trans. Signal Process., vol. 57, pp.
2479–2493, Jul. 2009.
[37] H. W. Kuhn, “The Hungarian method for the assignment problem,”
Naval research logistics quarterly, vol. 2, pp. 83–97, 1955.
[38] M. J. D. Powell, “On search directions for minimization algorithms,”
Math. Programming, vol. 4, pp. 193–201, Dec. 1973.
[39] L. Grippo and M. Sciandrone, “On the convergence of the block
nonlinear Gauss–Seidel method under convex constraints,” Operations
research letters, vol. 26, pp. 127–136, Apr. 2000.

Thomas Lundgaard Hansen received the B.Sc. and
M.Sc. (cum laude) in electrical engineering from
Aalborg University, Denmark in 2011 and 2014,
respectively. Since 2014 he has been pursuing the
Ph.D. degree in wireless communication at Aalborg
University. During 2013 and 2015 he was a visiting
scholar at the University of California, San Diego,
USA. He is the recipient of the best student paper
award (1st place) at the 2014 IEEE Sensor Array and
Multichannel Signal Processing workshop and also
received an award from IDA Efondet for his master’s
thesis. His research interests include signal processing, machine learning,
optimization and wireless communication.

Bernard Henri Fleury (M’97–SM’99) received the
Diplomas in Electrical Engineering and in Mathematics in 1978 and 1990 respectively and the
Ph.D. Degree in Electrical Engineering in 1990 from
the Swiss Federal Institute of Technology Zurich
(ETHZ), Switzerland. Since 1997, he has been with
the Department of Electronic Systems, Aalborg University, Denmark, as a Professor of Communication
Theory. From 2000 till 2014 he was Head of Section,
first of the Digital Signal Processing Section and
later of the Navigation and Communications Section.
From 2006 to 2009, he was partly affiliated as a Key Researcher with the
Telecommunications Research Center Vienna (ftw.), Austria. Prof. Fleury’s
research interests cover numerous aspects within communication theory, signal
processing, and machine learning, mainly for wireless communication systems
and networks.

Bhaskar D. Rao (S’80–M’83–SM’91–F’00) is currently a Distinguished Professor in the Electrical and
Computer Engineering department and the holder
of the Ericsson endowed chair in Wireless Access
Networks at the University of California, San Diego.
Prof. Rao was elected fellow of IEEE in 2000 and
is the recipient of the 2016 IEEE Signal Processing
Society Technical Achievement Award. Prof. Rao’s
interests are in the areas of digital signal processing,
estimation theory, and optimization theory, with applications to digital communications, speech signal
processing, and biomedical signal processing.

