# Digital Beamforming Robust to Time-Varying Carrier Frequency Offset

Shuang Li, Payam Nayeri, and Michael B. Wakin SL is with the Department of
Mathematics, University of California, Los Angeles, CA 90095. Email:
<EMAIL_ADDRESS>PN and MBW are with the Department of Electrical
Engineering, Colorado School of Mines, Golden, CO 80401. Email:
<EMAIL_ADDRESS>
###### Abstract

Adaptive interference cancellation is rapidly becoming a necessity for our
modern wireless communication systems, due to the proliferation of wireless
devices that interfere with each other. To cancel interference, digital
beamforming algorithms adaptively adjust the weight vector of the antenna
array, and in turn its radiation pattern, to minimize interference while
maximizing the desired signal power. While these algorithms are effective in
ideal scenarios, they are sensitive to signal corruptions. In this work, we
consider the case when the transmitter and receiver in a communication system
cannot be synchronized, resulting in a carrier frequency offset that corrupts
the signal. We present novel beamforming algorithms that are robust to signal
corruptions arising from this time-variant carrier frequency offset. In
particular, we bring in the Discrete Prolate Spheroidal Sequences (DPSS’s) and
propose two atomic-norm-minimization (ANM)-based methods in both 1D and 2D
frameworks to design a weight vector that can be used to cancel interference
when there exist unknown time-varying frequency drift in the pilot and
interferer signals. Both algorithms do not assume a pilot signal is known.
Noting that solving ANM optimization problems via semi-definite programs can
be a computational burden, we also present a novel fast algorithm to
approximately solve our 1D ANM optimization problem. Finally, we confirm the
benefits of our proposed algorithms and show the advantages over existing
approaches with a series of experiments.

###### Index Terms:

Array processing, atomic norm minimization (ANM), carrier frequency offset,
digital beamforming, discrete prolate spheroidal sequences (DPSS’s),
interference cancellation

##  I Introduction

Modern communication systems need to operate in a crowded spectrum where
multiple transceivers communicate in close physical proximity and in many
cases across very close channels. Interference has become a major issue,
degrading communication reliability and throughput [1], [2]. To survive this
hostile environment, communication systems need to adapt to the dynamic
environment by adaptively cancelling interference sources and maximizing
signal-to-interference-plus-noise ratio (SINR) [3], [4]. A digital beamformer
can adaptively adjust its weight vector using an array processing algorithm,
which in turn shapes its radiation pattern and can place pattern nulls in the
direction of interference sources, while maximizing the desired signal power
[5]. Digital beamforming (DBF) is thus seen as a necessity for future wireless
systems. While multiple approaches to cost-effective DBF hardware have been
proposed over the years [6, 7, 8], one of the major challenges remaining is
creating robust DBF algorithms that can operate effectively in the hostile
wireless signal environment.

A wide variety of algorithms have been developed for DBF [4, 9, 10], each
differing in what is assumed known and how signals are modeled. Very broadly,
these algorithms can be categorized into statistical versus deterministic
methods. Statistical methods [4, 9] utilize a stochastic model for the desired
signal, interferers, and/or noise, while deterministic methods [10] do not
rely on stochastic models and tend to involve eigendecompositions and related
concepts from linear algebra. DBF methods can also be categorized into blind
versus non-blind methods. Although there is no universal definition of blind
methods [11, 12, 10, 13], the gist is that such methods do not rely on
explicitly knowing aspects such as the desired signal, its spectral
properties, the array geometry, the signal direction, and so on. Depending on
what is assumed, blind beamforming problems can also have connections to blind
source separation problems [10, 13]. Non-blind methods may require explicit
knowledge of a pilot signal or its spectral properties as well as other
knowledge of the array geometry.

In this paper, we consider the DBF problem in the following context. Our goal
is to design a weight vector (see Section III) that allows us to localize and
separate a desired signal (“the signal” for short) from unwanted interferers
impinging on an array. Our treatment can be characterized as:

  * •

Deterministic: We do not assume a stochastic model for the signal,
interferers, or noise, nor do we require long-term observations for estimating
auto- or cross-correlation statistics.

  * •

Narrowband: We assume that the signal and interferers have narrowband spectra
relative to the sampling frequency. Sinusoids, which are a natural choice for
pilot signals (discussed below) are one possibility for such narrowband
signals, but sinusoidal structure is not explicitly required.

  * •

Blind: Our algorithms do not require a known pilot signal. We present two DBF
algorithms in this paper. Our 1D algorithm does not require the array antennas
to be equally spaced, nor does it require the antenna positions to be known.
Our 2D algorithm does require the antenna positions to be known and equally
spaced, but it can handle the case where the signal and interferers have
overlapping temporal frequencies.

Our approach to solving this DBF problem is inspired by several existing
approaches, which we describe next before giving further details on our
approach.

The minimum variance distortionless response (MVDR) beamformer is a
prototypical non-blind statistical beamformer [4]. Modeling the signal and
noise as stochastic processes, and given knowledge of the correlation
statistics (spectral matrices) for these processes in addition to knowledge of
the signal direction, MVDR is designed to output the highest possible SINR. As
such it is considered the benchmark for comparing the performance of
beamforming algorithms [4, 14]. In practice, however, some or all of these
quantities may not be available. As we discuss in Section II-A, when a known
pilot signal and the empirical data covariance are used to estimate the signal
direction and spectral properties, MVDR reduces to a technique known as sample
matrix inversion (SMI).111Other algorithms in this vein include least mean
square (LMS) and recursive least square (RLS). Both converge to the optimal
solution but are generally slower than SMI [4, 15]. Computing the SMI
beamformer requires only the empirical data matrix and knowledge of the pilot
signal (but not the array geometry or pilot signal direction). As such, it can
also be interpreted as a semi-blind deterministic method [11].

DBF algorithms such as SMI can effectively cancel interferers as long as the
knowledge of the incoming pilot signal is sufficiently accurate. However, in
non-ideal scenarios, signals can become corrupted and cause these algorithms
fail. To address these concerns, DBF algorithms need to be robust to causes
that corrupt the signal. Again, the literature on robust DBF methods is vast.
Several algorithms have been developed to address beamformer mismatches in
particular in the form of steering vector (SV) errors. SV errors arise from
direction-of-arrival mismatch, or perturbed arrays with errors in element
location, phase, or gain, and cause signal corruptions that degrades the
system performance. The most widely adopted solution to these problems is
diagonal loading and equivalent approaches [14, 16, 17, 18, 19]. Diagonal
loading provides robustness to SV errors by effectively designing for a higher
white noise level than is actually present, and as such is usually practical
when signal-to-interference ratio (SIR) is more important than SINR. To
mitigate some of these issues, variable loading has been introduced; it can
improve robustness to SV errors while maintaining a desired SIR or SINR [18,
19].

A different category of beamformer mismatch, which is the focus of this work,
occurs when the transmitter and receiver cannot be synchronized and thus a
carrier frequency offset exists between the anticipated pilot signal and that
which arrives at the array. In [15], it was shown that when it is not possible
to measure and correct the carrier frequency offset, algorithms such as SMI
will fail. A solution to this problem using genetic algorithms was
demonstrated in [15]. More recently, a new mathematical approach based on
atomic norm minimization (ANM) was proposed for this problem [20, 21]. The ANM
approach is naturally capable of extracting the offset frequencies that cause
signal corruption, and as such, robustness is achieved with maximum SINR and
without any degradation to other system performance characteristics. Both of
these algorithms provided a robust solution to beamforming when transmitter
and receiver are unsynchronized, which is a realistic non-ideal scenario
occurring more frequently in our ever growing modern wireless infrastructure.
This robust beamforming problem is categorically different than what has been
studied in the past, and to the best of our knowledge, no solution to this
problem is available outside [15, 20, 21].

While the works [15, 20, 21] provided a robust solution to the carrier
frequency offset problem in unsynchronized systems, they only considered the
time-invariant case. In reality, however, frequency offset is a dynamic
problem, and for practical implementation a robust time-variant solution is
needed. In this work, we accommodate an unknown time-varying frequency offset
in the pilot signal. Indeed, we drop entirely the requirement that the pilot
signal be known, and we assume only that the pilot signal is narrowband. As we
show (see Section II-C), this assumption accommodates (but does not require)
the common choice of a sinusoidal pilot signal, and it allows for a moderate
amount of unknown time-varying drift in the carrier frequency of this
sinusoid. In this way, our treatment is more general (blind) and robust
compared to classical SMI.

Our approach for solving the blind beamforming problem relies on an
optimization framework using ANM. ANM allows for a signal or data matrix to be
decomposed into sparse components arising from a certain atomic set of
prototype signals or matrices. In compressive sensing, finite atomic sets of
on-grid sinusoids have been used for recovering signals from partial
information [22, 23]. In line spectrum estimation, continuous atomic sets of
off-grid sinusoids have been used for precisely identifying the active
frequencies in a finite-length sampled signal without suffering from the
“leakage” effects that confront conventional Fourier-based methods [24]. In
array processing, such atomic sets have been used to provide joint sparse
frequency decompositions of multi-channel data matrices [25, 26].
Unfortunately, the atomic sets in these works are not rich enough to represent
the signals encountered when time-varying frequency drift is present.

In this paper, we expand the array processing ANM frameworks using a novel
atomic set for data matrices. Our work expands on the joint sparse frequency
model [25, 26] by modulating each sinusoid by an unknown vector in a low-
dimensional basis of Discrete Prolate Spheroidal Sequences (DPSS’s) [27, 28,
29]. DPSS’s have been shown to provide excellent compressed representations of
sampled narrowband signals. We therefore bring DPSS’s into the ANM framework
to capture the unknown time-varying frequency drift component in the pilot and
interferer signals.

We present two novel ANM frameworks in this paper. Our first “1D” framework
does not require the antennas to be equally spaced, and moreover it does not
require the antenna positions to be known. Our second “2D” framework does
require the antenna positions to be known and equally spaced, but it can
handle the case where the pilot signal and interferers have overlapping
temporal frequencies. The 2D framework exploits sparse structure in both the
temporal and angular frequency domains, whereas the 1D framework exploits
sparse sinusoidal structure only in the temporal frequency domain. Both
frameworks can accommodate time-varying frequency drift component in the pilot
and interferer signals.

Finally, we note that most ANM optimization problems can be formulated as a
semi-definite program (SDP), solving this SDP can be a computational burden.
In this paper, we also present a novel fast algorithm for approximately
solving our 1D ANM optimization problem. We provide extensive experimental
tests to confirm the benefits of our 1D and 2D ANM frameworks compared to
conventional SMI, and we illustrate the computational speedup provided by our
fast 1D algorithm. We leave the work of developing a fast algorithm for the 2D
ANM method as our future work.

Notations. Throughout this paper, we use superscripts ⊤, ∗, H, and † to denote
transpose, conjugate, conjugate transpose, and pseudo-inverse, respectively.
We use $\odot$ and $\otimes$ to denote the elementwise multiplication and
Kronecker product, respectively. For a positive integer $K$, denote set
$[K]\triangleq\\{1,\cdots,K\\}$. We use $\mathbf{1}_{L}\in\mathbb{R}^{L}$ and
$\mathbf{I}_{L}\in\mathbb{R}^{L\times L}$ to represent a length-$L$ vector
with entries being all 1’s and an $L\times L$ identity matrix, respectively.
Finally, we let
$\mathbf{P}_{\mathbf{Q}}\triangleq\mathbf{I}-\mathbf{Q}\mathbf{Q}^{H}$ denote
the orthogonal projection onto the orthogonal complement of the column space
of an orthogonal matrix $\mathbf{Q}$.

The remainder of this paper is organized as follows. We introduce some
background on SMI, ANM, and DPSS in Section II and formulate the problem of
interest in Section III. We then present the proposed 1D method
“ANM+DPSS+SMI”, fast 1D method “IVDST+DPSS+SMI”, and 2D method
“2D-ANM+DPSS+SMI” in Sections IV, V, and VI, respectively. We conduct a series
of experiments to test the proposed methods in Section VII. Finally, we
conclude this work and discuss a few possible future directions in Section
VIII.

##  II Preliminaries

###  II-A Sample matrix inversion (SMI)

Consider an $N$-element linear antenna array. The optimal MVDR beamformer
computes the complex valued weighting vector $\bm{w}\in\mathbb{C}^{N}$ for the
array elements such that the interference signals are minimized subject to a
constraint [4]. The algorithm minimizes
$\mathbb{E}\left\\{\bm{w}^{H}\mathbf{S}_{x}\bm{w}\right\\}$ subject to
$\bm{w}^{H}\textbf{asv}=1$, where the constraint ensures that a plane wave
with a known array steering vector, asv, will not be distorted. Here,
$\mathbb{E}$ is the expected value, and $\mathbf{S}_{x}$ is the signal
spectral matrix. In practice $\mathbf{S}_{x}$ is not available but the maximum
likelihood of the matrix can be estimated from incoming data as
$\mathbf{S}_{x}=\frac{1}{M}\sum_{m=1}^{M}\bm{x}_{m}\bm{x}_{m}^{H}$, where
$\\{\bm{x}_{m}\\}_{m=1}^{M}$ denote the signal snapshots and $M$ is the number
of snapshots. The array steering vector is
$\textbf{asv}=e^{jk_{0}\sin(\theta_{s})\bm{q}}$, where
$\bm{q}\in\mathbb{R}^{1\times N}$ denotes the element positions, $\theta_{s}$
is the signal direction, and $k_{0}$ is the plane wave wavenumber. In SMI,
where a pilot signal is available, the signal steering vector, $\bm{v}_{s}$,
is used in place of the array steering vector, and its maximum likelihood is
estimated by the cross-correlation between the signal, $\bm{x}$, and the known
pilot signal, $\bm{s}\in\mathbb{C}^{M}$, as
$\bm{v}_{s}=\frac{1}{M}\sum_{m=1}^{M}\bm{x}_{m}s_{m}^{H},$ where $s_{m}$ is
the $m$-th entry of $\bm{s}$. The SMI element weights are then given by
$\bm{w}^{H}\triangleq\frac{\bm{v}_{s}^{H}\mathbf{S}_{x}^{-1}}{\bm{v}_{s}^{H}\mathbf{S}_{x}^{-1}\bm{v}_{s}}$.
In this equation, the presence of the $\mathbf{S}_{x}^{-1}$ term gives rise to
the name sample matrix inversion. By rearranging the signal snapshots
$\\{\bm{x}_{m}\\}_{m=1}^{M}$ as columns of a data matrix $\mathbf{X}^{H}$, it
follows that the SMI element weights can be rewritten as

$\bm{w}=(\mathbf{X}^{H}\mathbf{X})^{-1}\mathbf{X}^{H}\bm{s}=\mathbf{X}^{\dagger}\bm{s},$
(1)

where we ignored the normalization factor
$\bm{v}_{s}^{H}\mathbf{S}_{x}^{-1}\bm{v}_{s}$. Though we arrived at the
solution (1) using statistical arguments, we note that
$\bm{w}=\mathbf{X}^{\dagger}\bm{s}$ can also be interpreted as the least-
squares solution to the purely deterministic problem of finding the weight
vector that best solves $\mathbf{X}\bm{w}=\bm{s}$, i.e., best localizes the
known pilot signal $\bm{s}$ [11].

###  II-B Atomic norm minimization (ANM)

The atomic norm is a generalization of the $\ell_{1}$ norm (commonly used in
compressive sensing and sparse signal recovery) and nuclear norm (commonly
used in low-rank matrix recovery) from finite, discrete dictionaries to
infinite, continuously-parameterized dictionaries [30, 24]. As in $\ell_{1}$
norm minimization and nuclear norm minimization, atomic norm minimization
(ANM) can promote sparsity when solving an inverse problem. Thus, it has been
widely studied in many signal processing applications, including line spectral
estimation and array signal processing, for localizing the off-grid temporal
or angular frequencies of sinusoidal components [30, 31, 24, 25, 26, 32, 33,
34].

To be more precise, consider the following spectrally sparse signal
$\bm{x}\in\mathbb{C}^{M}$ with $K$ different active frequencies

$\displaystyle\bm{x}=\sum_{k=1}^{K}c_{k}\bm{a}(f_{k})\in\mathbb{C}^{M},$ (2)

where the $c_{k}$’s denote complex coefficients, and

$\displaystyle\bm{a}(f)\triangleq[1\leavevmode\nobreak\ e^{j2\pi
f}\leavevmode\nobreak\ \cdots\leavevmode\nobreak\ e^{j2\pi
f(M-1)}]^{\top}\in\mathbb{C}^{M}$ (3)

is a vector containing $M$ uniform samples of a complex exponential signal
with frequency $f\in[0,1)$. One can then define

$\mathcal{A}\triangleq\\{\bm{a}(f):f\in[0,1)\\}$

as an atomic set that consists of all possible such complex exponentials. Note
that the spectrally sparse signal $\bm{x}$ in (2) can be represented with only
$K$ atoms from the above atomic set $\mathcal{A}$. The corresponding induced
atomic norm is then defined as

$\|\bm{x}\|_{\mathcal{A}}\triangleq\inf\left\\{\sum_{k=1}^{K}|c_{k}|:\bm{x}=\sum_{k=1}^{K}c_{k}\bm{a}(f_{k})\right\\}.$

###  II-C Discrete prolate spheroidal sequences (DPSS’s)

For integers $M$ and $L\in\\{1,2,\dots,M\\}$ and a positive number
$W<\frac{1}{2}$, let $\mathbf{S}_{M,W}\in\mathbb{R}^{M\times L}$ denote a
matrix containing the first $L$ of the $M$-dimensional discrete prolate
spheroidal sequence (DPSS) vectors with digital half bandwidth $W$ [27, 28,
29]. DPSS’s are distinctive in that each column of $\mathbf{S}_{M,W}$, if zero
padded, has a Discrete-Time Fourier Transform (DTFT) highly concentrated in
the normalized frequency domain $[-W,W]$. In fact, the first column vector is
precisely the time-limited signal (with $M$ nonzero entries) whose DTFT is
optimally concentrated in $[-W,W]$, the second column vector is the time-
limited signal whose DTFT is most concentrated in $[-W,W]$ subject to being
orthogonal to the first column, etc.

DPSS’s are particularly valuable for providing efficient representations of
sampled sinusoids and narrowband signals. For example, as indicated in the
following theorem, any complex exponential signal $\bm{a}(f)$ defined in (3)
with $f\in[-W,W]$ can be well represented in the column span of the DPSS basis
$\mathbf{S}_{M,W}$ when $L\approx 2MW$ [28, 29].

######  Theorem II.1.

[29] Fix $W\in(0,\frac{1}{2})$ and $\epsilon\in(0,1)$. Let
$L=2MW(1+\epsilon)$. Then there exist constants $C_{1},C_{2}$ such that

$\begin{split}||\mathbf{P}_{\mathbf{S}_{M,W}}\bm{a}(f)||_{2}\leq
C_{1}M^{5/4}e^{-C_{2}M},\leavevmode\nobreak\ \forall\leavevmode\nobreak\
f\in[-W,W]\end{split}$

holds for all $M\geq M_{0}$. Here, $C_{1},C_{2}$ are numerical constants that
may depend on $W$ and $\epsilon$.

In summary, this result states that ($i$) the effective dimensionality of the
subspace spanned by $\\{\bm{a}(f)\\}_{f\in[-W,W]}$ is $2MW$, and ($ii$) the
first $L\approx 2MW$ DPSS vectors provide a basis for this subspace. Sinusoids
concentrated in a narrow band away from baseband, e.g., $f\in[a,b]$, can be
equally well-captured by taking $W=\frac{b-a}{2}$ to be the half bandwidth and
modulating each of the DPSS vectors by $\bm{a}(f^{o})$ where
$f^{0}=\frac{a+b}{2}$ is the center frequency. Finally, the high quality of
DPSS approximations applies not only to sampled sinusoids, but also to sampled
narrowband signals (since these can be viewed as linear combinations of
sampled sinusoids within a narrow range of frequencies). As discussed more
formally in [28, 29], the finite-length sample vectors arising from sampling
random baseband (or bandpass) analog signals can be well approximated using a
DPSS (or modulated DPSS) basis. Again the effective dimensionality is
$L\approx 2MW$, where $M$ is the length of the sample vector and $W$ is the
digital half bandwidth.

##  III Problem Formulation

In this work, we consider a conventional linear array with $N$ antenna
elements. We collect $M$ snapshots on each antenna element and model the data
matrix as

$\displaystyle\mathbf{X}^{\star}=\sum_{k=1}^{K}\bm{s}_{k}\textbf{asv}(\theta_{k})\in\mathbb{C}^{M\times
N},$ (4)

where $\bm{s}_{1}\in\mathbb{C}^{M}$ denotes the desired signal and
$\\{\bm{s}_{k}\in\mathbb{C}^{M}\\}_{k=2}^{K}$ denote the interferers. Here,

$\displaystyle\textbf{asv}(\theta_{k})\triangleq
e^{jk_{0}\sin(\theta_{k})\bm{q}}$

is a length-$N$ row vector representing the array steering vector associated
with the desired signal (if $k=1$) or interferers (if $k>1$). $\theta_{1}$ and
$\\{\theta_{k}\\}_{k=2}^{K}$ denote the angle to the desired signal and
interferers, respectively, $\bm{q}\in\mathbb{R}^{1\times N}$ denotes the
element positions, and $k_{0}$ is the plane wave wavenumber.

We do not assume the desired signal, the interferers, or any of their
directions are known. Moreover, our 1D algorithm does not require the array
antennas to be equally spaced, nor does it require the element positions
$\bm{q}$ to be known.

Given the data matrix $\mathbf{X}^{\star}$, our goal is to construct a weight
vector $\bm{w}$ such that $\mathbf{X}^{\star}\bm{w}\approx\bm{s}_{1}$ as
closely as possible. Inspired by SMI (which assumes $\bm{s}_{1}$ is known),
our approach involves first estimating $\bm{s}_{1}$ from $\mathbf{X}^{\star}$
and then plugging this estimate (which we call $\widetilde{\bm{s}}_{1}$) into
the SMI expression:

$\displaystyle\bm{w}={\mathbf{X}^{\star}}^{\dagger}\widetilde{\bm{s}}_{1}.$
(5)

To facilitate the estimation of $\bm{s}_{1}$ from $\mathbf{X}^{\star}$, we
assume that the desired signal and all of the interferers are narrowband
signals whose DTFT’s are each concentrated in a small subset of the digital
bandwidth $[-\frac{1}{2},\frac{1}{2}]$. More specifically, employing the
DPSS’s discussed in Section II-C, we assume that each $\bm{s}_{k}$ can be
modeled as follows:

$\bm{s}_{k}\approx\bm{a}(f_{k}^{o})\odot(\mathbf{S}_{M,W}\boldsymbol{\alpha}_{k})\leavevmode\nobreak\
\text{for}\leavevmode\nobreak\ k=1,2,\dots,K,$ (6)

where $f_{k}^{o}$ is a possibly unknown carrier frequency, $\bm{a}(f_{k}^{o})$
is a sampled complex exponential with frequency $f_{k}^{o}$ as defined in (3),
$\mathbf{S}_{M,W}\in\mathbb{R}^{M\times L}$ is an $L$-dimensional DPSS basis
(we discuss the choice of $W$ and $L\approx 2MW$), and
$\boldsymbol{\alpha}_{k}\in\mathbb{C}^{L}$ denote the corresponding unknown
coefficients. This narrowband assumption is quite flexible and accommodates
the following cases:

  * •

Case 1: $\bm{s}_{k}$ is a pure sinusoid (complex exponential) with a fixed
frequency. In this case, (6) holds with equality by taking $f_{k}^{o}$ to be
the frequency of the sinusoid and letting $\boldsymbol{\alpha}_{k}=0$.

  * •

Case 2: $\bm{s}_{k}$ is a sinusoid (such as a pilot signal when $k=1$),
modulated by a frequency offset that drifts slightly over time. In this case,
$\bm{s}_{k}$ can be viewed as a frequency modulated (FM) signal. Carson’s
bandwidth rule for FM signals [35] guarantees that a majority of the signal
power will concentrate in a narrow band around the signal’s carrier frequency,
with the half bandwidth $W$ proportional to the peak frequency deviation plus
the temporal bandwidth of the modulating signal. Therefore,
$\bm{a}(f_{k}^{o})$ will again capture the dominant sinusoidal component, and
following the discussion in Section II-C,
$\mathbf{S}_{M,W}\boldsymbol{\alpha}_{k}$ will approximately capture the
deviations from this frequency.

  * •

Case 3: $\bm{s}_{k}$ is a bandpass signal, oversampled with respect to its
Nyquist rate. In this case, the DTFT of $\bm{s}_{k}$ will concentrate in some
subset $[a,b]$ of the digital bandwidth $[-\frac{1}{2},\frac{1}{2}]$.
Following the discussion in Section II-C, the approximation in (6) then holds
by taking $f_{k}^{o}=\frac{a+b}{2}$ to be the midpoint of this band and $W$ to
be greater than or equal to its half width $\frac{b-a}{2}$.

  * •

Case 4: $\bm{s}_{k}$ is a bandpass signal, oversampled with respect to its
Nyquist rate, modulated by a frequency offset that drifts slightly over time.
The approximation in (6) then holds by combining the arguments from Case 2 and
Case 3: because multiplication in time corresponds to convolution in
frequency, the product of two narrowband signals will remain (relatively)
narrowband.

The reader will notice that we do not focus on quantifying the necessary
degree of accuracy in the approximation (6) in this work. Nevertheless, this
assumption is quite broad. We use numerical experiments in Section VII to
confirm that this approximation is sufficiently valid to enable the success of
our proposed technique in a range of interesting cases.

Finally, we note that the above idea of using (5) to design a weight vector
has already been studied in [20, 21] under the assumption that each
$\bm{s}_{k}$ is a pure sinusoid modulated by a fixed (not time-varying)
frequency offset. Those works propose two ANM-based methods to estimate the
frequency offset from either $\mathbf{X}^{\star}$ or its noisy measurements.
However, as is shown in the simulations, the weight vector obtained by
directly applying this idea cannot cancel the interferers when there exist
time-varying frequency offsets. To address this problem, we extend the ANM
formulation using the novel assumption in (6).

##  IV Proposed Method: “ANM+DPSS+SMI”

We now describe our procedure for constructing the estimate
$\widetilde{\bm{s}}_{1}$ based on the data matrix $\mathbf{X}^{\star}$. We
will estimate $\bm{s}_{1}$ by estimating both $f_{1}^{o}$ and
$\boldsymbol{\alpha}_{1}$ and plugging these estimates into (6). Indeed, we
can estimate any and all of the frequencies $f_{k}^{o}$ and vectors
$\boldsymbol{\alpha}_{k}$ using the procedure below; we assume only that,
among these estimates, the correct one can be associated with the desired
signal (the $k=1$ term can be correctly selected).

Inspired by a recent work on multi-band signal recovery [36], we rewrite the
data matrix in (4) as

$\displaystyle\mathbf{X}^{\star}$
$\displaystyle=\sum_{k=1}^{K}\bm{s}_{k}\textbf{asv}(\theta_{k})\overset{\text{\char
172}}{\approx}\sum_{k=1}^{K}\left(\bm{a}(f_{k}^{o})\odot(\mathbf{S}_{M,W}\boldsymbol{\alpha}_{k})\right)\textbf{asv}(\theta_{k})$
$\displaystyle=\sum_{k=1}^{K}\|\boldsymbol{\alpha}_{k}\|_{2}\|\textbf{asv}(\theta_{k})\|_{2}((\mathbf{S}_{M,W}\text{sign}(\boldsymbol{\alpha}_{k}))$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
\odot\bm{a}(f_{k}^{o}))\text{sign}(\textbf{asv}(\theta_{k}))$
$\displaystyle\overset{\text{\char
173}}{=}\sum_{k=1}^{K}c_{k}\left((\mathbf{S}_{M,W}\widetilde{\boldsymbol{\alpha}}_{k})\odot\bm{a}(f_{k}^{o})\right)\bm{b}_{k}^{H}$
$\displaystyle=\sum_{k=1}^{K}c_{k}[((\mathbf{S}_{M,W}\widetilde{\boldsymbol{\alpha}}_{k})\odot\bm{a}(f_{k}^{o}))b_{1k}^{*}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\cdots((\mathbf{S}_{M,W}\widetilde{\boldsymbol{\alpha}}_{k})\odot\bm{a}(f_{k}^{o}))b_{Nk}^{*}]$
$\displaystyle\overset{\text{\char
174}}{=}\sum_{k=1}^{K}c_{k}\left[\left(\mathbf{S}_{M,W}\odot\left(\bm{a}(f_{k}^{o})\widetilde{\boldsymbol{\alpha}}_{k}^{\top}b_{1k}^{*}\right)\right)\mathbf{1}_{L}\right.$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\cdots\left.\left(\mathbf{S}_{M,W}\odot\left(\bm{a}(f_{k}^{o})\widetilde{\boldsymbol{\alpha}}_{k}^{\top}b_{Nk}^{*}\right)\right)\mathbf{1}_{L}\right]$
$\displaystyle=\left[\left(\mathbf{S}_{M,W}\odot\left(\sum_{k=1}^{K}c_{k}\bm{a}(f_{k}^{o})\widetilde{\boldsymbol{\alpha}}_{k}^{\top}b_{1k}^{*}\right)\right)\mathbf{1}_{L}\right.$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\left.\cdots\left(\mathbf{S}_{M,W}\odot\left(\sum_{k=1}^{K}c_{k}\bm{a}(f_{k}^{o})\widetilde{\boldsymbol{\alpha}}_{k}^{\top}b_{Nk}^{*}\right)\right)\mathbf{1}_{L}\right],$

where the sign vector is defined as
$\text{sign}(\bm{x})\triangleq\frac{\bm{x}}{\|\bm{x}\|_{2}}$. ① follows from
our assumption (6). ② follows by denoting
$c_{k}=\|\boldsymbol{\alpha}_{k}\|_{2}\|\textbf{asv}(\theta_{k})\|_{2}$,
$\widetilde{\boldsymbol{\alpha}}_{k}=\text{sign}(\boldsymbol{\alpha}_{k})$,
and $\bm{b}_{k}^{H}=\text{sign}(\textbf{asv}(\theta_{k}))$. Note that
$\\{\widetilde{\boldsymbol{\alpha}}_{k}\\}_{k=1}^{K}$ and
$\\{\bm{b}_{k}\\}_{k=1}^{K}$ are normalized vectors with
$\|\widetilde{\boldsymbol{\alpha}}_{k}\|_{2}=\|\bm{b}_{k}\|_{2}=1$. We use
$b_{nk}$ to denote the $n$-th entry of $\bm{b}_{k}$.
$\mathbf{1}_{L}\in\mathbb{R}^{L}$ is a length-$L$ vector with entries being
all 1’s. Finally, ③ follows from the following equalities

$\displaystyle(\mathbf{S}_{M,W}\widetilde{\boldsymbol{\alpha}}_{k})\odot\bm{a}(f_{k}^{o})$
$\displaystyle=\text{diag}(\bm{a}(f_{k}^{o}))\mathbf{S}_{M,W}\widetilde{\boldsymbol{\alpha}}_{k}$
$\displaystyle=\text{diag}(\bm{a}(f_{k}^{o}))\mathbf{S}_{M,W}\text{diag}(\widetilde{\boldsymbol{\alpha}}_{k})\mathbf{1}_{L}$
$\displaystyle=\left(\mathbf{S}_{M,W}\odot\left(\bm{a}(f_{k}^{o})\widetilde{\boldsymbol{\alpha}}_{k}^{\top}\right)\right)\mathbf{1}_{L}.$

The above observation inspires us to define a tensor
$\mathcal{X}^{\star}\in\mathbb{C}^{M\times L\times N}$ with the $n$-th frontal
slice being
$\mathcal{X}^{\star}_{::n}=\sum_{k=1}^{K}c_{k}\bm{a}(f_{k}^{o})\widetilde{\boldsymbol{\alpha}}_{k}^{\top}b_{nk}^{*}$,
namely,

$\displaystyle\mathcal{X}^{\star}\triangleq\sum_{k=1}^{K}c_{k}\mathbf{A}(f_{k}^{o})\circledast\mathbf{H}_{k}$
(7)

with
$\mathbf{A}(f_{k}^{o})\triangleq[\bm{a}(f_{k}^{o})\cdots\bm{a}(f_{k}^{o})]\in\mathbb{C}^{M\times
N}$ and
$\mathbf{H}_{k}\triangleq\widetilde{\boldsymbol{\alpha}}_{k}\bm{b}_{k}^{H}\in\mathbb{C}^{L\times
N}$. Here, we use $\mathbf{A}\circledast\mathbf{B}\in\mathbb{C}^{M\times
L\times N}$ to denote the “reshaped Khatri-Rao product”, which is defined as
$[\mathbf{A}\circledast\mathbf{B}]_{::n}\triangleq\bm{a}_{n}\bm{b}_{n}^{\top}$
for any two matrices
$\mathbf{A}=[\bm{a}_{1}\cdots\bm{a}_{N}]\in\mathbb{C}^{M\times N}$ and
$\mathbf{B}=[\bm{b}_{1}\cdots\bm{b}_{N}]\in\mathbb{C}^{L\times N}$. By
defining a linear operator $\mathcal{L}:\mathbb{C}^{M\times L\times
N}\rightarrow\mathbb{C}^{M\times N}$ as

$\displaystyle\mathcal{L}(\mathcal{X})\triangleq[(\mathbf{S}_{M,W}\odot\mathcal{X}_{::1})\mathbf{1}_{L}\cdots(\mathbf{S}_{M,W}\odot\mathcal{X}_{::N})\mathbf{1}_{L}],$
(8)

we can rewrite the data matrix $\mathbf{X}^{\star}$ as the linear measurements
of a third-order tensor $\mathcal{X}^{\star}$ obtained from the above linear
operator $\mathcal{L}:\mathbb{C}^{M\times L\times
N}\rightarrow\mathbb{C}^{M\times N}$, namely,

$\displaystyle\mathbf{X}^{\star}=\mathcal{L}(\mathcal{X}^{\star}).$

We now define an atomic set as

$\displaystyle\mathcal{A}\triangleq\left\\{\mathbf{A}(f^{o})\circledast\mathbf{H}:\leavevmode\nobreak\
f^{o}\in[0,1),\leavevmode\nobreak\ \mathbf{H}\in\mathbb{C}^{L\times
N},\leavevmode\nobreak\ \|\mathbf{H}\|_{F}=1\right\\}.$

The induced atomic norm is then defined as

$\displaystyle\|\mathcal{X}\|_{\mathcal{A}}\triangleq\inf\left\\{\sum_{k=1}^{K}c_{k}:\mathcal{X}=\sum_{k=1}^{K}c_{k}\mathbf{A}(f_{k}^{o})\circledast\mathbf{H}_{k},c_{k}\geq
0\right\\},$

which is approximately equivalent to the optimal value of the following
semidefinite program (SDP) [36]:

$\displaystyle\inf_{\begin{subarray}{c}\mathbf{U}\in\mathbb{C}^{M\times N}\\\
\mathcal{W}\in\mathbb{C}^{L\times L\times N}\end{subarray}}$
$\displaystyle\frac{1}{2}u+\frac{1}{2}\sum_{n=1}^{N}\text{Tr}(\mathcal{W}_{::n})$
$\displaystyle\operatorname*{s.t.}\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ $
$\displaystyle\left[\begin{array}[]{cc}\text{Toep}(\bm{u}_{n})&\mathcal{X}_{::n}\\\
\mathcal{X}_{::n}^{H}&\mathcal{W}_{::n}\end{array}\right]\succeq\mathbf{0},\leavevmode\nobreak\
u_{1n}=u,\leavevmode\nobreak\ \forall\leavevmode\nobreak\ n\in[N].$

Here, we use Tr$(\cdot)$ to denote the trace of a square matrix and
Toep$(\bm{u})$ to denote the Hermitian Toeplitz matrix with the vector
$\bm{u}$ as its first column. $\bm{u}_{n}$ and $u_{1n}$ denote the $n$-th
column of $\mathbf{U}$ and the first entry in $\bm{u}_{n}$, respectively.

It can be seen that the tensor defined in (7) can be represented with only $K$
atoms from the above atomic set. To promote sparsity and recover the tensor
$\mathcal{X}^{\star}$ from its linear measurements
$\mathbf{X}^{\star}=\mathcal{L}(\mathcal{X}^{\star})$, we propose to solve the
following ANM problem

$\displaystyle\min_{\mathcal{X}}\leavevmode\nobreak\
\|\mathcal{X}\|_{\mathcal{A}}\quad\operatorname*{s.t.}\leavevmode\nobreak\
\|\mathbf{X}^{\star}-\mathcal{L}(\mathcal{X})\|_{F}\leq\varepsilon,$ (9)

where $\varepsilon$ is a parameter that controls the data fidelity term. Note
that the above ANM (9) is approximately equivalent to the following SDP

$\displaystyle\inf_{\begin{subarray}{c}\mathbf{U}\in\mathbb{C}^{M\times N}\\\
\mathcal{W}\in\mathbb{C}^{L\times L\times N}\\\
\mathcal{X}\in\mathbb{C}^{M\times L\times N}\end{subarray}}$
$\displaystyle\frac{1}{2}u+\frac{1}{2}\sum_{n=1}^{N}\text{Tr}(\mathcal{W}_{::n})$
(10) $\displaystyle\operatorname*{s.t.}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak$
$\displaystyle\left[\begin{array}[]{cc}\text{Toep}(\bm{u}_{n})&\mathcal{X}_{::n}\\\
\mathcal{X}_{::n}^{H}&\mathcal{W}_{::n}\end{array}\right]\succeq\mathbf{0},\leavevmode\nobreak\
u_{1n}=u,\leavevmode\nobreak\ \forall\leavevmode\nobreak\ n\in[N],$
$\displaystyle\|\mathbf{X}^{\star}-\mathcal{L}(\mathcal{X})\|_{F}\leq\varepsilon,$

which can be solved by any off-the-shelf SDP solver, such as CVX [37]. As
primal-dual algorithms are used in CVX, it can return both the primal solution
and the dual solution.

Denote $\mathcal{Q}^{\star}$ as the dual solution to (9). We construct a dual
polynomial as

$\displaystyle
q(f^{o})\triangleq\|\boldsymbol{\Psi}(f^{o},\mathcal{Q}^{\star})\|_{F}$ (11)

with

$\displaystyle\boldsymbol{\Psi}(f^{o},\mathcal{Q}^{\star})=\left[\begin{array}[]{ccc}{\mathcal{Q}^{\star}_{::1}}^{H}\bm{a}(f^{o})&\cdots&{\mathcal{Q}^{\star}_{::N}}^{H}\bm{a}(f^{o})\end{array}\right]\in\mathbb{C}^{L\times
N}.$

As is shown in [36], one can identify $\\{f_{k}^{o}\\}_{k=1}^{K}$ by
localizing the places where $q(f^{o})=1$. However, in this work, as we will
see in the simulations, the dual polynomial has several plateaus. Therefore,
we will use the k-means method to cluster the $f^{o}$ with
$q(f^{o})\geq\gamma_{0}$ and get the estimated frequencies
$\widetilde{f}_{k}^{o}$.

Denote $\widehat{\mathcal{X}}$ as the primal optimal solution to (9). Note
that

$\displaystyle\widehat{\mathcal{X}}_{::1}$
$\displaystyle=\sum_{k=1}^{K}c_{k}\bm{a}(\widetilde{f}_{k}^{o})\widehat{\widetilde{\boldsymbol{\alpha}}}_{k}^{\top}\widehat{b}_{1k}^{*}$
$\displaystyle=[\bm{a}(\widetilde{f}_{1}^{o})\cdots\bm{a}(\widetilde{f}_{K}^{o})]\left[\begin{array}[]{ccc}c_{1}\widehat{b}_{11}^{*}&&\\\
&\ddots&\\\
&&c_{K}\widehat{b}_{1K}^{*}\end{array}\right]\left[\begin{array}[]{c}\widehat{\widetilde{\boldsymbol{\alpha}}}_{1}^{\top}\\\
\vdots\\\
\widehat{\widetilde{\boldsymbol{\alpha}}}_{K}^{\top}\end{array}\right].$

Denote
$\mathbf{A}_{f}\triangleq[\bm{a}(\widetilde{f}_{1}^{o})\cdots\bm{a}(\widetilde{f}_{K}^{o})]$.
Once we have the estimated frequencies
$\\{\widetilde{f}_{k}^{o}\\}_{k=1}^{K}$, we can then estimate the sign of the
coefficient $\widehat{\widetilde{\boldsymbol{\alpha}}}_{1}$ as

$\displaystyle\text{sign}(\widehat{\widetilde{\boldsymbol{\alpha}}}_{1})=\frac{(\mathbf{A}_{f}^{\dagger}\widehat{\mathcal{X}}_{::1})^{\top}\bm{e}_{1}}{\|(\mathbf{A}_{f}^{\dagger}\widehat{\mathcal{X}}_{::1})^{\top}\bm{e}_{1}\|_{2}},$

where $\bm{e}_{1}\in\mathbb{R}^{K}$ is the first column of the $K\times K$
identity matrix $\mathbf{I}_{K}$. Recall that
$\widetilde{\boldsymbol{\alpha}}_{k}=\text{sign}(\boldsymbol{\alpha}_{k})$.
Therefore, we get

$\displaystyle\text{sign}(\widehat{\boldsymbol{\alpha}}_{1})=\frac{(\mathbf{A}_{f}^{\dagger}\widehat{\mathcal{X}}_{::1})^{\top}\bm{e}_{1}}{\|(\mathbf{A}_{f}^{\dagger}\widehat{\mathcal{X}}_{::1})^{\top}\bm{e}_{1}\|_{2}}.$
(12)

Finally, we estimate $\bm{s}_{1}$ as

$\displaystyle\widetilde{\bm{s}}_{1}=\bm{a}(\widetilde{f}_{1}^{o})\odot(\mathbf{S}_{M,W}\text{sign}(\widehat{\boldsymbol{\alpha}}_{1}))$
(13)

and we plug this estimate into (5) to estimate the weight vector $\bm{w}$. We
summarize the proposed “ANM+DPSS+SMI” method in Algorithm 1.

Algorithm 1 ANM+DPSS+SMI

1:procedure Input(the data matrix $\mathbf{X}^{\star}$ and the number of DPSS
basis vectors $L$)

2: create DPSS basis $\mathbf{S}_{M,W}\in\mathbb{R}^{M\times L}$

3: compute the primal solution $\widehat{\mathcal{X}}$ and the dual solution
$\mathcal{Q}^{\star}$ by solving the SDP (10)

4: form the dual polynomial $q(f^{o})$ as in (11)

5: use the k-means method to cluster the $f^{o}$ with $q(f^{o})\geq\gamma_{0}$
and get $\widetilde{f}_{k}^{o}$; set $k=1$ according to the desired frequency

6: compute $\text{sign}(\widehat{\boldsymbol{\alpha}}_{1})$ according to (12)

7: estimate $\bm{s}_{1}$ as in (13) and construct the weight vector as in (5)

8: return the weight vector $\bm{w}$

9:end procedure

######  Remark IV.1.

(Computational complexity.) Note that most SDP solvers, including SDPT3 [38],
use the interior-point method and need to solve a system of linear equations
when computing the Newton direction, which can be very expensive for large-
size problems. Here, the $N$ positive semi-definite (PSD) constraints in (10)
can be combined into one PSD constraint on a big block diagonal matrix of size
$N(M+L)\times N(M+L)$. Therefore, the overall computational complexity of
solving (10) with SDP solvers is $\mathcal{O}(N^{3.5}(M+L)^{3.5})$, which
limits the use of the proposed “ANM+DPSS+SMI” method in large-scale problems.

In Remark IV.1, one may notice that the combined PSD constraint is on a block
diagonal matrix, rather than a dense matrix. Thus, a natural question arises:

  * Is it possible that the computational complexity of solving (10) is less than $\mathcal{O}(N^{3.5}(M+L)^{3.5})$ due to the sparsity that exists in the PSD constraint?

To the best of our knowledge, there is no closed-form expression for the
complexity. However, we are able to empirically answer this question by
conducting the following experiment. In particular, we collect the running
time needed for solving the following SDP

$\min\leavevmode\nobreak\ x\leavevmode\nobreak\ \leavevmode\nobreak\
\operatorname*{s.t.}\leavevmode\nobreak\ \leavevmode\nobreak\
\mathbf{A}_{x}\succeq\mathbf{0},$ (14)

where $\mathbf{A}_{x}\in\mathbb{R}^{N_{s}\times N_{s}}$ is a matrix with
entries from the set $\\{0,x,\sqrt{N_{s}}\\}$. First, we test the case when
$\mathbf{A}_{x}$ is a dense matrix, i.e.,

$\displaystyle\mathbf{A}_{x}=\left[\begin{array}[]{cccc}x&\sqrt{N_{s}}&\cdots&\sqrt{N_{s}}\\\
\sqrt{N_{s}}&x&\cdots&\sqrt{N_{s}}\\\ \vdots&\vdots&\ddots&\vdots\\\
\sqrt{N_{s}}&\sqrt{N_{s}}&\cdots&x\end{array}\right].$ (19)

The time needed for solving (14) with $\mathbf{A}_{x}$ given in (19) for a
variety of $N_{s}$ is illustrated in Figure 1 (the blue line with circles).
Second, we test several cases when $\mathbf{A}_{x}$ is a block diagonal
matrix, specifically with $N_{b}$ of blocks on the diagonal where
$N_{b}=10,\leavevmode\nobreak\ 20,\leavevmode\nobreak\ 50$ and $100$. We
generate each block on the diagonal as an
$\frac{N_{s}}{N_{b}}\times\frac{N_{s}}{N_{b}}$ matrix with diagonal entries
being $x$ and off-diagonal entries being $\sqrt{N_{s}}$. The running times are
presented as red lines in Figure 1. It can be seen that solving the SDP (14)
with a sparse PSD constraint matrix is only slightly faster when using a dense
PSD constraint matrix. The running time for both cases is approximately on the
same order as the constraint matrix size, i.e., $\mathcal{O}(N_{s}^{2.5})$,
which is better than the complexity described in Remark IV.1.

Figure 1: The time (in seconds) needed for solving SDP (14) with varying size
of PSD constraint matrix (denoted as $N_{s}$).

We conclude that the answer to the above question is “no”, namely, that
sparsity in the PSD constraint matrix will not significantly reduce the
computational complexity of solving an SDP problem. To reduce the
computational complexity, we introduce a fast version of the “ANM+DPSS+SMI”
method in next section.

##  V Proposed Method: “IVDST+DPSS+SMI”

Recall that a key step in the proposed “ANM+DPSS+SMI” method is to compute the
dual solution. Though CVX can return the dual solution by solving the primal
SDP (10), one can also directly solve the dual problem of (10) to get the dual
solution. With some elementary calculations, we obtain the dual problem of
(10) as follows222To simplify the derivation, we replace the inequality
constraint
$\|\mathbf{X}^{\star}-\mathcal{L}(\mathcal{X})\|_{F}\leq\varepsilon$ with an
equality constraint $\mathbf{X}^{\star}=\mathcal{L}(\mathcal{X})$, i.e., we
set $\varepsilon=0$.

$\displaystyle\sup_{\begin{subarray}{c}\mathcal{Q}\in\mathbb{C}^{M\times
L\times N}\\\ \mathcal{H}\in\mathbb{C}^{M\times M\times N}\end{subarray}}$
$\displaystyle\langle\mathbf{X}^{\star},\mathcal{L}(\mathcal{Q})\rangle_{\mathbb{R}}$
(20) $\displaystyle\operatorname*{s.t.}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak$
$\displaystyle\left[\begin{array}[]{cc}\mathcal{H}_{::n}&-\mathcal{Q}_{::n}\\\
-\mathcal{Q}_{::n}^{H}&\mathbf{I}_{L}\end{array}\right]\succeq\mathbf{0},$
$\displaystyle\sum_{m=1}^{M-j}\mathcal{H}_{m(m+j)n}=0,\leavevmode\nobreak\
\forall\leavevmode\nobreak\ j\in[M-1],\leavevmode\nobreak\
\forall\leavevmode\nobreak\ n\in[N],$
$\displaystyle\sum_{n=1}^{N}\sum_{m=1}^{M}\mathcal{H}_{mmn}=1.$

Note that the above dual problem (20) is also an SDP. Therefore, it can be
solved by any off-the-shelf SDP solver just as the primal SDP (10). However,
as is discussed in Remark IV.1, solving SDP problems directly with SDP solvers
can result in high computational complexity. To reduce the computational
complexity, the authors in [39] propose a fast method based on Alternating
Direction Method of Multipliers (ADMM) [40, 41] to solve the SDP. In [42], the
authors reformulate the SDP as a conic program with much fewer dual variables
to further reduce the computational cost. A fast iterative Vandermonde
decomposition and shrinkage-thresholding (IVDST) algorithm based on the
accelerated proximal gradient (APG) technique is developed in [43] to solve
the SDP. Recently, this algorithm is also extended to the two-dimensional (2D)
framework in a problem of 2D grid-free compressive beamforming [44].

Inspired by the IVDST algorithm, we propose a faster method, named
“IVDST+DPSS+SMI”, to speed up the “ANM+DPSS+SMI” method introduced in Section
IV. Rather than solving the primal SDP (10), we focus on its dual problem (20)
in this section. Note that the dual problem (20) can be rewritten as the
following unconstrained problem

$\displaystyle\min_{\mathcal{Q},\mathcal{H}}\leavevmode\nobreak\
h(\mathcal{Q})+g(\mathcal{Q},\mathcal{H})$

with

$\displaystyle h(\mathcal{Q})$
$\displaystyle\triangleq-\langle\mathbf{X}^{\star},\mathcal{L}(\mathcal{Q})\rangle_{\mathbb{R}},$
$\displaystyle g(\mathcal{Q},\mathcal{H})$ $\displaystyle\triangleq
g_{1}(\mathcal{Q},\mathcal{H})+g_{2}(\mathcal{H})+g_{3}(\mathcal{H}),$
$\displaystyle g_{1}(\mathcal{Q},\mathcal{H})$
$\displaystyle\triangleq\lambda_{1}\sum_{n=1}^{N}\mathcal{C}(\mathbf{Z}(\mathcal{Q}_{::n},\mathcal{H}_{::n})),$
$\displaystyle g_{2}(\mathcal{H})$
$\displaystyle\triangleq\lambda_{2}\left(\sum_{n=1}^{N}\sum_{m=1}^{M}\mathcal{H}_{mmn}-1\right)^{2},$
$\displaystyle g_{3}(\mathcal{H})$
$\displaystyle\triangleq\lambda_{3}\sum_{n=1}^{N}\sum_{j=1}^{M-1}\left(\sum_{m=1}^{M-j}\mathcal{H}_{m(m+j)n}\right)^{2},$

where $\\{\lambda_{i}\\}_{i=1}^{3}$ are regularization parameters.
$\mathbf{Z}(\mathcal{Q}_{::n},\mathcal{H}_{::n})$ is a matrix defined as

$\displaystyle\mathbf{Z}(\mathcal{Q}_{::n},\mathcal{H}_{::n})\triangleq\left[\begin{array}[]{cc}\mathcal{H}_{::n}&-\mathcal{Q}_{::n}\\\
-\mathcal{Q}_{::n}^{H}&\mathbf{I}_{L}\end{array}\right].$

$\mathcal{C}(\mathbf{Z})$ is an indicator function used to enforce the PSD
constraint and is defined as

$\displaystyle\mathcal{C}(\mathbf{Z})\triangleq\begin{cases}0,\quad&\text{if}\leavevmode\nobreak\
\mathbf{Z}\leavevmode\nobreak\ \text{is PSD,}\\\
\infty,&\text{otherwise.}\end{cases}$

Next, we introduce how to use the idea of APG to estimate the dual solution
$\mathcal{Q}^{\star}$, which is then used to design a weight vector $\bm{w}$.

(1) Initialization:

According to the Schur complement condition,
$\left[\begin{array}[]{cc}\mathcal{H}_{::n}&-\mathcal{Q}_{::n}\\\
-\mathcal{Q}_{::n}^{H}&\mathbf{I}_{L}\end{array}\right]\succeq\mathbf{0}$ and
$\mathbf{I}_{L}\succ\mathbf{0}$ imply that

$\displaystyle\mathcal{H}_{::n}-\mathcal{Q}_{::n}\mathcal{Q}_{::n}^{H}\succeq\mathbf{0}.$

Therefore, we initialize
$\mathcal{Q}_{1}=\mathcal{Q}_{0}\in\mathbb{C}^{M\times L\times N}$ as a
complex Gaussian random tensor, and initialize
$\mathcal{H}_{1}=\mathcal{H}_{0}\in\mathbb{C}^{M\times M\times N}$ as

$\displaystyle\\{\mathcal{H}_{1}\\}_{::n}=\\{\mathcal{H}_{0}\\}_{::n}=\\{\mathcal{Q}_{0}\\}_{::n}\\{\mathcal{Q}_{0}\\}_{::n}^{H},\quad\forall\leavevmode\nobreak\
n\in[N].$ (21)

(2) Smoothing:

In each iteration, we add a momentum term to accelerate the convergence of the
gradient vector. In particular, the updates in the $i$-th iteration are given
as

$\displaystyle\overline{\mathcal{Q}}_{i}=\mathcal{Q}_{i}+\frac{t_{i-1}-1}{t_{i}}(\mathcal{Q}_{i}-\mathcal{Q}_{i-1}),$
(22)
$\displaystyle\overline{\mathcal{H}}_{i}=\mathcal{H}_{i}+\frac{t_{i-1}-1}{t_{i}}(\mathcal{H}_{i}-\mathcal{H}_{i-1}),$

where $t_{i}=\frac{1+\sqrt{4t_{i-1}^{2}+1}}{2}$ with $t_{0}=1$.

(3) Gradient descent:

After smoothing, we update the parameters $\mathcal{Q}$ along the gradient
descent direction of $h(\mathcal{Q})$ with an appropriate small stepsize
$\eta$. $\mathcal{H}$ is fixed in this step since $h(\mathcal{Q})$ only
depends on $\mathcal{Q}$. Then, the updates in the $i$-th iteration are given
as

$\displaystyle\mathcal{Q}_{g}$
$\displaystyle=\overline{\mathcal{Q}}_{i}-\eta\nabla
h(\mathcal{Q})=\overline{\mathcal{Q}}_{i}+\eta\mathcal{L}^{*}(\mathbf{X}^{\star}),$
(23) $\displaystyle\mathcal{H}_{g}$
$\displaystyle=\overline{\mathcal{H}}_{i},$

where $\mathcal{L}^{*}:\mathbb{C}^{M\times N}\rightarrow\mathbb{C}^{M\times
L\times N}$ denotes the adjoint operator of the linear operator $\mathcal{L}$
defined in (8). In particular, we have

$\displaystyle[\mathcal{L}^{*}(\mathbf{X}^{\star})]_{::n}=\mathbf{S}_{M,W}\odot(\bm{x}^{\star}_{n}\mathbf{1}_{L}^{H}),\quad\forall\leavevmode\nobreak\
n\in[N],$ (24)

where $\bm{x}^{\star}_{n}$ denotes the $n$-th column of $\mathbf{X}^{\star}$.

(4) Proximal mapping:

In this step, we use the following proximal operator

$\displaystyle(\mathcal{Q}_{i+1},\mathcal{H}_{i+1})=\text{prox}_{\eta
g}(\mathcal{Q}_{g},\mathcal{H}_{g})$ $\displaystyle=$
$\displaystyle\arg\min_{\mathcal{Q},\mathcal{H}}\left\\{g(\mathcal{Q},\mathcal{H})+\frac{1}{2\eta}\|(\mathcal{Q},\mathcal{H})-(\mathcal{Q}_{g},\mathcal{H}_{g})\|_{F}^{2}\right\\}.$

Note that there is no analytical solution to the above optimization problem,
therefore, we use alternating projection method to approximate its solution.
In particular, for any $n\in[N]$, we first update $\mathcal{H}$ with

$\displaystyle\widetilde{\mathcal{H}}_{mmn}$
$\displaystyle=\frac{\\{\mathcal{H}_{g}\\}_{mmn}}{\sum_{n=1}^{N}\sum_{m=1}^{M}\\{\mathcal{H}_{g}\\}_{mmn}},$
(25) $\displaystyle\widetilde{\mathcal{H}}_{m(m+j)n}$
$\displaystyle=\\{\mathcal{H}_{g}\\}_{m(m+j)n}-\frac{1}{M-j}\sum_{m=1}^{M-j}\\{\mathcal{H}_{g}\\}_{m(m+j)n},$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\forall\leavevmode\nobreak\ j\in[M-1],$

which ensures

$\displaystyle\sum_{n=1}^{N}\sum_{m=1}^{M}\widetilde{\mathcal{H}}_{mmn}=1,\text{
and }$
$\displaystyle\sum_{m=1}^{M-j}\widetilde{\mathcal{H}}_{m(m+j)n}=0,\leavevmode\nobreak\
\forall\leavevmode\nobreak\ j\in[M-1],\leavevmode\nobreak\
\forall\leavevmode\nobreak\ n\in[N],$

as is required in the constraints of (20).

Next, to promote the PSD constraint, we evaluate

$\displaystyle\mathbf{Z}_{n}=\left[\begin{array}[]{cc}\widetilde{\mathcal{H}}_{::n}&-\\{\mathcal{Q}_{g}\\}_{::n}\\\
-\\{\mathcal{Q}_{g}\\}_{::n}^{H}&\mathbf{I}_{L}\end{array}\right],\quad\forall\leavevmode\nobreak\
n\in[N]$ (28)

and compute its eigen-decomposition

$\displaystyle(\mathbf{V}_{n},\boldsymbol{\Sigma}_{n})=\text{eig}(\mathbf{Z}_{n}),\quad\forall\leavevmode\nobreak\
n\in[N].$

To make $\\{\mathbf{Z}_{n}\\}_{n=1}^{N}$ PSD, we only keep the positive
eigenvalues and the corresponding eigenvectors. Then, we get

$\displaystyle\widetilde{\mathbf{Z}}_{n}=\mathbf{V}_{n}\max(\boldsymbol{\Sigma}_{n},0)\mathbf{V}_{n}^{H},\leavevmode\nobreak\
\forall\leavevmode\nobreak\ n\in[N].$ (29)

With the obtained PSD matrices $\\{\widetilde{\mathbf{Z}}_{n}\\}_{n=1}^{N}$,
we update $\mathcal{H}$ and $\mathcal{Q}$ with

$\displaystyle\\{\mathcal{H}_{i+1}\\}_{::n}$
$\displaystyle=\widetilde{\mathbf{Z}}_{n}(1:M,1:M),$ (30)
$\displaystyle\\{\mathcal{Q}_{i+1}\\}_{::n}$
$\displaystyle=-\widetilde{\mathbf{Z}}_{n}(1:M,M+1:M+L),\quad\forall\leavevmode\nobreak\
n\in[N].$

(5) Computing weight:

Denoting $I$ as the number of iterations, then $\mathcal{Q}_{I}$ is an
estimation of the dual solution $\mathcal{Q}^{\star}$. Given
$\mathcal{Q}_{I}$, we can then construct a dual polynomial as in (11) and use
the k-means method to cluster the $f^{o}$ with $q(f^{o})\geq\gamma_{0}$ and
get the estimated frequencies $\widetilde{f}_{k}^{o}$. As is shown in (12), to
compute $\text{sign}(\widehat{\boldsymbol{\alpha}}_{1})$, one needs to first
estimate $\widehat{\mathcal{X}}_{::1}$, which is the primal solution to the
SDP (10). However, since we focus on solving the dual SDP (20), we do not have
the primal solution. Fortunately, we can use
$[\mathcal{L}^{*}(\mathbf{X}^{\star})]_{::1}$ as a surrogate to
$\widehat{\mathcal{X}}_{::1}$. Then, we can compute
$\text{sign}(\widehat{\boldsymbol{\alpha}}_{1})$ as

$\displaystyle\text{sign}(\widehat{\boldsymbol{\alpha}}_{1})=\frac{(\mathbf{A}_{f}^{\dagger}[\mathcal{L}^{*}(\mathbf{X}^{\star})]_{::1})^{\top}\bm{e}_{1}}{\|(\mathbf{A}_{f}^{\dagger}[\mathcal{L}^{*}(\mathbf{X}^{\star})]_{::1})^{\top}\bm{e}_{1}\|_{2}}$
(31)

with
$\mathbf{A}_{f}=[\bm{a}(\widetilde{f}_{1}^{o})\cdots\bm{a}(\widetilde{f}_{K}^{o})]$.
Finally, $\widetilde{\bm{s}}_{1}$ can be computed as in (13) and the weight
vector $\bm{w}$ can be computed using (5).

We summarize the proposed “IVDST+DPSS+SMI” method in Algorithm 2.

Algorithm 2 IVDST+DPSS+SMI

1:procedure Input(the data matrix $\mathbf{X}^{\star}$, the number of DPSS
basis vectors $L$, the stepsize $\eta$, and the number of iterations $I$)

2: create DPSS basis $\mathbf{S}_{M,W}\in\mathbb{R}^{M\times L}$

3: compute $\mathcal{L}^{*}(\mathbf{X}^{\star})$ according to (24)

4: initialize $\mathcal{Q}_{1}=\mathcal{Q}_{0}$ as a complex Gaussian random
tensor, and initialize $\mathcal{H}_{1}=\mathcal{H}_{0}$ as in (21)

5: for $i=1,...,I$ do

6: smoothing: compute $\overline{\mathcal{Q}}_{i}$ and
$\overline{\mathcal{H}}_{i}$ via (22)

7: gradient descent: compute $\mathcal{Q}_{g}$ and $\mathcal{H}_{g}$ via (23)

8: proximal mapping: compute $\widetilde{\mathcal{H}}$ via (25), evaluate
$\mathbf{Z}_{n}$ via (28) and compute its eigen-composition to get
$(\mathbf{V}_{n},\boldsymbol{\Sigma}_{n})$, then compute
$\widetilde{\mathbf{Z}}_{n}$ via (29), finally, update $\mathcal{Q}_{i+1}$ and
$\mathcal{H}_{i+1}$ via (30)

9: end for

10: return dual solution $\mathcal{Q}_{I}$

11: given $\mathcal{Q}_{I}$, form the dual polynomial $q(f^{o})$ as in (11)

12: use the k-means method to cluster the $f^{o}$ with
$q(f^{o})\geq\gamma_{0}$ and get $\widetilde{f}_{k}^{o}$; set $k=1$ according
to the desired frequency

13: compute $\text{sign}(\widehat{\boldsymbol{\alpha}}_{1})$ according to (31)

14: estimate $\bm{s}_{1}$ as in (13) and construct the weight vector as in (5)

15: return the weight vector $\bm{w}$

16:end procedure

######  Remark V.1.

(Computational complexity.) Compared with the proposed “ANM+DPSS+SMI” method,
the most expensive step in the “IVDST+DPSS+SMI” method is computing the eigen-
decomposition of $\mathbf{Z}_{n},\leavevmode\nobreak\
\forall\leavevmode\nobreak\ n\in[N]$. Therefore, the overall complexity is on
the order of $\mathcal{O}(N(M+L)^{2})$.

##  VI Extension to 2D Case

In this section, we take the structure of $\textbf{asv}(\theta)$ into
consideration and define a two dimensional (2D) atomic set as

$\displaystyle\mathcal{A}_{2D}\triangleq\left\\{\mathbf{A}(f^{o})\circledast(\widetilde{\boldsymbol{\alpha}}\textbf{asv}(\theta)):\leavevmode\nobreak\
f^{o}\in[0,1),\leavevmode\nobreak\
\widetilde{\boldsymbol{\alpha}}\in\mathbb{C}^{L},\right.$ (32)
$\displaystyle\left.\|\widetilde{\boldsymbol{\alpha}}\|_{2}=1,\leavevmode\nobreak\
\theta\in[-90^{\circ},90^{\circ}]\right\\}.$

The induced 2D atomic norm is then defined as

$\displaystyle\|\mathcal{X}\|_{\mathcal{A}_{2D}}\triangleq\inf\left\\{\sum_{k=1}^{K}c_{k}:\right.\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak$ (33) $\displaystyle\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\left.\mathcal{X}=\sum_{k=1}^{K}c_{k}\mathbf{A}(f_{k}^{o})\circledast(\widetilde{\boldsymbol{\alpha}}_{k}\textbf{asv}(\theta_{k})),c_{k}\geq
0\right\\},$

As in Section IV, one can solve the following 2D-ANM problem

$\displaystyle\min_{\mathcal{X}}\leavevmode\nobreak\
\|\mathcal{X}\|_{\mathcal{A}_{2D}}\quad\operatorname*{s.t.}\leavevmode\nobreak\
\|\mathbf{X}^{\star}-\mathcal{L}(\mathcal{X})\|_{F}\leq\varepsilon$ (34)

to construct the weight vector.

Denote
$\langle\mathcal{Q},\mathcal{X}\rangle_{\mathbb{R}}=\text{Re}\\{\langle\mathcal{Q},\mathcal{X}\rangle\\}$
as the real inner product between two tensors. The dual norm of
$\|\mathcal{X}\|_{\mathcal{A}_{2D}}$ is then defined as

$\displaystyle\|\mathcal{Q}\|_{\mathcal{A}_{2D}}^{*}=$
$\displaystyle\sup_{\|\mathcal{X}\|_{\mathcal{A}_{2D}}\leq
1}\langle\mathcal{Q},\mathcal{X}\rangle_{\mathbb{R}}$ $\displaystyle=$
$\displaystyle\sup_{\begin{subarray}{c}f^{o}\in[0,1),\|\widetilde{\boldsymbol{\alpha}}\|_{2}=1\\\
\theta\in[-90^{\circ},90^{\circ}]\end{subarray}}\langle\mathcal{Q},\mathbf{A}(f^{o})\circledast(\widetilde{\boldsymbol{\alpha}}\textbf{asv}(\theta))\rangle_{\mathbb{R}}$
$\displaystyle=$
$\displaystyle\sup_{\begin{subarray}{c}f^{o}\in[0,1),\|\widetilde{\boldsymbol{\alpha}}\|_{2}=1\\\
\theta\in[-90^{\circ},90^{\circ}]\end{subarray}}\sum_{n=1}^{N}\langle\mathcal{Q}_{::n},\bm{a}(f^{o})e^{jk_{0}\sin(\theta)\bm{q}_{n}}\widetilde{\boldsymbol{\alpha}}^{\top}\rangle_{\mathbb{R}}$
$\displaystyle=$
$\displaystyle\sup_{\begin{subarray}{c}f^{o}\in[0,1),\|\widetilde{\boldsymbol{\alpha}}\|_{2}=1\\\
\theta\in[-90^{\circ},90^{\circ}]\end{subarray}}\left\langle\widetilde{\boldsymbol{\alpha}}^{*},\sum_{n=1}^{N}\mathcal{Q}_{::n}^{H}\bm{a}(f^{o})e^{jk_{0}\sin(\theta)\bm{q}_{n}}\right\rangle_{\mathbb{R}}$
$\displaystyle=$ $\displaystyle\sup_{\begin{subarray}{c}f^{o}\in[0,1)\\\
\theta\in[-90^{\circ},90^{\circ}]\end{subarray}}\left\|\sum_{n=1}^{N}\mathcal{Q}_{::n}^{H}\bm{a}(f^{o})e^{jk_{0}\sin(\theta)\bm{q}_{n}}\right\|_{2}$
$\displaystyle=$ $\displaystyle\sup_{\begin{subarray}{c}f^{o}\in[0,1)\\\
\theta\in[-90^{\circ},90^{\circ}]\end{subarray}}q_{2D}(f^{o},\theta),$

where

$\displaystyle
q_{2D}(f^{o},\theta)\triangleq\left\|\sum_{n=1}^{N}\mathcal{Q}_{::n}^{H}\bm{a}(f^{o})e^{jk_{0}\sin(\theta)\bm{q}_{n}}\right\|_{2}$
(35)

is defined as the dual polynomial and $\bm{q}_{n}$ denotes the $n$-th element
position.

Next, we show that the above 2D atomic norm defined in (33) can be
approximated by the optimal value of the following semidefinite program (SDP):

$\displaystyle\inf_{\begin{subarray}{c}\mathcal{T}\in\mathbb{C}^{(2M-1)\times(2L-1)\times
N}\\\ \bm{t}\in\mathbb{C}^{N}\end{subarray}}$
$\displaystyle\frac{1}{2MN}\sum_{n=1}^{N}\text{Tr}(\SS(\mathcal{T}_{::n}))+\frac{1}{2N}\sum_{n=1}^{N}t_{n}$
(36) $\displaystyle\operatorname*{s.t.}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak$
$\displaystyle\left[\begin{array}[]{cc}\SS(\mathcal{T}_{::n})&\bm{x}_{n}\\\
\bm{x}_{n}^{H}&t_{n}\end{array}\right]\succeq\mathbf{0},$
$\displaystyle\bm{x}_{n}=\text{vec}(\mathcal{X}_{::n}^{\top}),\leavevmode\nobreak\
\forall\leavevmode\nobreak\ n\in[N],$

where $\SS(\mathbf{T})\in\mathbb{C}^{ML\times ML}$ is a block Toeplitz matrix
generated by $\mathbf{T}\in\mathbb{C}^{(2M-1)\times(2L-1)}$. To be more
precise, denote $T_{m,l}$ with $-M<m<M$ and $-L<l<L$ as the $(m,l)$-th entry
of $\mathbf{T}$. Denote $\mathbf{T}_{m}\in\mathbb{C}^{L\times L}$ as a
Toeplitz matrix constructed from the $m$-th row of $\mathbf{T}$, i.e.,

$\displaystyle\mathbf{T}_{m}=\left[\begin{array}[]{cccc}T_{m,0}&T_{m,-1}&\cdots&T_{m,-(L-1)}\\\
T_{m,1}&T_{m,0}&\cdots&T_{m,-(L-2)}\\\ \vdots&\vdots&\ddots&\vdots\\\
T_{m,L-1}&T_{m,L-2}&\cdots&T_{m,0}\\\ \end{array}\right].$

Then, the block Toeplitz matrix $\SS(\mathbf{T})$ is defined as

$\displaystyle\SS(\mathbf{T})=\left[\begin{array}[]{cccc}\mathbf{T}_{0}&\mathbf{T}_{-1}&\cdots&\mathbf{T}_{-(M-1)}\\\
\mathbf{T}_{1}&\mathbf{T}_{0}&\cdots&\mathbf{T}_{-(M-2)}\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathbf{T}_{M-1}&\mathbf{T}_{M-2}&\cdots&\mathbf{T}_{0}\\\
\end{array}\right].$

Denote SDP$(\mathcal{X})$ as the objective value corresponding to the optimal
solution of the above SDP (36). The following proposition shows that
SDP$(\mathcal{X})$ is a lower bound of $\|\mathcal{X}\|_{\mathcal{A}_{2D}}$.

######  Proposition 1.

With the 2D-atomic set and 2D-atomic norm defined in (32) and (33), we have

$\displaystyle\text{SDP}(\mathcal{X})\leq\|\mathcal{X}\|_{\mathcal{A}_{2D}}$

for any $\mathcal{X}\in\mathbb{C}^{M\times L\times N}$.

###### Proof.

Let

$\displaystyle\mathcal{X}=\sum_{k=1}^{K}c_{k}\mathbf{A}(f_{k}^{o})\circledast(\widetilde{\boldsymbol{\alpha}}_{k}\textbf{asv}(\theta_{k}))$

with $f^{o}_{k}\in[0,1),\leavevmode\nobreak\
\widetilde{\boldsymbol{\alpha}}_{k}\in\mathbb{C}^{L},\leavevmode\nobreak\
\|\widetilde{\boldsymbol{\alpha}}_{k}\|_{2}=1,\leavevmode\nobreak\
\theta_{k}\in[-90^{\circ},90^{\circ}]$. Then, for any $n\in[N]$, we have

$\displaystyle\mathcal{X}_{::n}=\sum_{k=1}^{K}c_{k}\bm{a}(f_{k}^{o})\widetilde{\boldsymbol{\alpha}}_{k}^{\top}\textbf{asv}_{n}(\theta_{k}),$

and

$\displaystyle\bm{x}_{n}=\text{vec}(\mathcal{X}_{::n}^{\top})=\sum_{k=1}^{K}c_{k}\bm{a}(f_{k}^{o})\otimes\widetilde{\boldsymbol{\alpha}}_{k}\textbf{asv}_{n}(\theta_{k}),$

where $\textbf{asv}_{n}(\theta_{k})$ denotes the $n$-th entry of
$\textbf{asv}(\theta_{k})$ and $\otimes$ denotes the Kronecker product.

For $t_{n}=\sum_{k=1}^{K}c_{k}$ and a block Toeplitz matrix

$\displaystyle\SS(\mathcal{T}_{::n})\\!=\\!\sum_{k=1}^{K}\\!c_{k}\\!\left[\bm{a}(f_{k}^{o})\\!\otimes\widetilde{\boldsymbol{\alpha}}_{k}\textbf{asv}_{n}(\theta_{k})\right]\left[\bm{a}(f_{k}^{o})\otimes\widetilde{\boldsymbol{\alpha}}_{k}\textbf{asv}_{n}(\theta_{k})\right]^{H}\\!\\!,$

we have

$\displaystyle\left[\begin{array}[]{cc}\SS(\mathcal{T}_{::n})&\bm{x}_{n}\\\
\bm{x}_{n}^{H}&t_{n}\end{array}\right]$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{K}c_{k}\left[\begin{array}[]{c}\bm{a}(f_{k}^{o})\otimes\widetilde{\boldsymbol{\alpha}}_{k}\textbf{asv}_{n}(\theta_{k})\\\
1\end{array}\right]\left[\begin{array}[]{c}\bm{a}(f_{k}^{o})\otimes\widetilde{\boldsymbol{\alpha}}_{k}\textbf{asv}_{n}(\theta_{k})\\\
1\end{array}\right]^{H}$ $\displaystyle\succeq$ $\displaystyle\mathbf{0}.$

Thus, the above $\SS(\mathcal{T}_{::n})$ and $t_{n}$ with any $n\in[N]$ are
feasible to the SDP (36). It follows that

$\displaystyle\text{SDP}(\mathcal{X})\leq$
$\displaystyle\frac{1}{2MN}\sum_{n=1}^{N}\text{Tr}(\SS(\mathcal{T}_{::n}))+\frac{1}{2N}\sum_{n=1}^{N}t_{n}$
$\displaystyle=$
$\displaystyle\frac{1}{2MN}\sum_{n=1}^{N}\sum_{k=1}^{K}c_{k}\|\bm{a}(f_{k}^{o})\otimes\widetilde{\boldsymbol{\alpha}}_{k}\textbf{asv}_{n}(\theta_{k})\|_{2}^{2}+\frac{1}{2}\sum_{k=1}^{K}c_{k}$
$\displaystyle=$
$\displaystyle\sum_{k=1}^{K}c_{k}=\|\mathcal{X}\|_{\mathcal{A}_{2D}}.$

∎

Then, we propose the following SDP to approximate the solution of 2D-ANM in
(34):

$\displaystyle\inf_{\begin{subarray}{c}\mathcal{T}\in\mathbb{C}^{(2M-1)\times(2L-1)\times
N}\\\ \bm{t}\in\mathbb{C}^{N}\\\ \mathcal{X}\in\mathbb{C}^{M\times L\times
N}\end{subarray}}$
$\displaystyle\frac{1}{2MN}\sum_{n=1}^{N}\text{Tr}(\SS(\mathcal{T}_{::n}))+\frac{1}{2N}\sum_{n=1}^{N}t_{n}$
(37) $\displaystyle\operatorname*{s.t.}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak$
$\displaystyle\left[\begin{array}[]{cc}\SS(\mathcal{T}_{::n})&\bm{x}_{n}\\\
\bm{x}_{n}^{H}&t_{n}\end{array}\right]\succeq\mathbf{0},$
$\displaystyle\bm{x}_{n}=\text{vec}(\mathcal{X}_{::n}^{\top}),\leavevmode\nobreak\
\forall\leavevmode\nobreak\ n\in[N],$
$\displaystyle\|\mathbf{X}^{\star}-\mathcal{L}(\mathcal{X})\|_{F}\leq\varepsilon.$

Again, solving the above SDP with CVX can return both the primal solution and
the dual solution. Given the dual solution, one can construct a dual
polynomial as in (35) and identify $\\{f_{k}^{o},\theta_{k}\\}_{k=1}^{K}$ by
localizing the places where $q_{2D}(f^{o},\theta)=1$, similar as in Section
IV, and finally construct a weight vector.

We summarize the proposed “2D-ANM+DPSS+SMI” method in Algorithm 3.

Algorithm 3 2D-ANM+DPSS+SMI

1:procedure Input(the data matrix $\mathbf{X}^{\star}$ and the number of DPSS
basis vectors $L$)

2: create DPSS basis $\mathbf{S}_{M,W}\in\mathbb{R}^{M\times L}$

3: compute the primal solution $\widehat{\mathcal{X}}$ and the dual solution
$\mathcal{Q}^{\star}$ by solving the SDP (37)

4: form the dual polynomial $q_{2D}(f^{o})$ as in (35)

5: use the k-means method to cluster the $f^{o}$ with
$q_{2D}(f^{o},\theta)\geq\gamma_{0}$ and get $\widetilde{f}_{k}^{o}$; set
$k=1$ according to the desired frequency

6: compute $\text{sign}(\widehat{\boldsymbol{\alpha}}_{1})$ according to (12)

7: estimate $\bm{s}_{1}$ as in (13) and construct the weight vector as in (5)

8: return the weight vector $\bm{w}$

9:end procedure

######  Remark VI.1.

(Advantages and disadvantages of “2D-ANM+DPSS+SMI”.) Unlike the 1D-based
methods, in the above proposed “2D-ANM+DPSS+SMI” method, one needs to know the
element positions $\bm{q}\in\mathbb{R}^{1\times N}$, and the element positions
need to be equispaced. However, the “2D-ANM+DPSS+SMI” method can successfully
identify the frequency components in the signals even when there exist two
identical frequencies from different directions while the 1D method must have
a certain separation between all frequencies. See Figure 9.

##  VII Numerical Simulations

In this section, we conduct a series of experiments to test the proposed three
methods and compare them with the “ANM+SMI” method proposed in [21] and the
classical SMI method [4]. We simulate a uniform linear array with $N=4$
elements and half-wavelength element spacing, namely, the element position
vector $\bm{q}=-(N-1)d/2:d:(N-1)d/2$ with $d=0.5$. We assume one desired
signal and two interferers; thus, $K=3$. The angles to the desired signal and
two interferers are set as $\theta_{1}=-20^{\circ}$, $\theta_{2}=-60^{\circ}$,
and $\theta_{3}=20^{\circ}$.

In all experiments, we generate each of $\bm{s}_{1},\bm{s}_{2},\bm{s}_{3}$ as
a sinusoid (complex exponential) modulated by a frequency offset that drifts
slightly over time. This corresponds to Case 2 in Section III and allows for
an illustrative and challenging set of experiments. The ground truth sinusoid
frequencies are given by $f_{1}^{o}=0.1$, $f_{2}^{o}=0.3$ and $f_{3}^{o}=0.5$.
We consider four different types of frequency offsets: (a) static frequency
offset, (b) linear frequency offset, (c) zigzag frequency offset, and (d)
random frequency offset, which are shown in Figure 2. The three colors in each
plot represent the time-varying offset relative to the ground truth
frequencies $f_{1}^{o}$, $f_{2}^{o}$, and $f_{3}^{o}$. We take $M=120$ uniform
time samples at each element and formulate the data matrix
$\mathbf{X}^{\star}$ as in (4).

For the proposed “IVDST+DPSS+SMI” method, we set the stepsize as $\eta=4$ and
decrease it by multiplying a factor of $0.99$ after every $50$ iterations. The
maximum number of iterations is set as $I=200$. We present the values of $L$
used to generate the DPSS basis in Table I.

The dual polynomials obtained from the “ANM+SMI” method, the “ANM+DPSS+SMI”
method, and the “IVDST+DPSS+SMI” method are shown in Figure 3. The
corresponding radiation pattern (array factor) is illustrated in Figure 4.
Figure 4 also includes the radiation pattern for “SMI” which corresponds to
the classical SMI method but accounting only for the sinusoidal component of
$\bm{s}_{1}$ and not its time-varying offset; that is, setting
$\bm{w}=(\mathbf{X}^{H}\mathbf{X})^{-1}\mathbf{X}^{H}\bm{s}=\mathbf{X}^{\dagger}\bm{a}(f_{1}^{o})$.
(Our methods do not assume explicit knowledge of $f_{1}^{o}$.)

From Figures 3 and 4, it can be seen that the “ANM+SMI” method [21] works very
well when the frequency offsets are static. However, our proposed methods
“IVDST+DPSS+SMI” and “ANM+DPSS+SMI” significantly outperform the “ANM+SMI” and
“SMI” methods when the frequency offsets are time-varying. Moreover, based on
the dual polynomials shown in Figure 3, it is much easier to cluster the
frequencies with the fast “IVDST+DPSS+SMI” method, as the dual polynomials
obtained from the “ANM+DPSS+SMI” method are sometimes relatively flat. To show
the efficiency of the proposed “IVDST+DPSS+SMI” method, we present the time
used by “IVDST+DPSS+SMI” and “ANM+DPSS+SMI” in Table II. As expected, the
“IVDST+DPSS+SMI” method runs much faster than the SDP-based “ANM+DPSS+SMI”
method.

(a) Static frequency offset

(b) Linear frequency offset

(c) Zigzag frequency offset

(d) Random frequency offset

Figure 2: Four types of frequency offsets used in the experiments. ($M=120$)

(a) Static frequency offset

(b) Linear frequency offset

(c) Zigzag frequency offset

(d) Random frequency offset

Figure 3: Dual polynomial obtained from the three methods in the four types of
frequency offsets. ($M=120$)

(a) Static frequency offset

(b) Linear frequency offset

(c) Zigzag frequency offset

(d) Random frequency offset

Figure 4: Interference cancellation with DBF in the four types of frequency offsets. ($M=120$) Table I: The value of $L$ used to generate the DPSS basis. | Static | Linear | Zigzag | Random
---|---|---|---|---
IVDST+DPSS+SMI | 7 | 2 | 2 | 13
ANM+DPSS+SMI | 7 | 10 | 10 | 13
2D-ANM+DPSS+SMI | 4 | 4 | 5 | 4
Table II: Running time (in seconds) of the “IVDST+DPSS+SMI” method and the “ANM+DPSS+SMI” method used in the first experiment. | Static | Linear | Zigzag | Random
---|---|---|---|---
IVDST+DPSS+SMI | 20.25 | 19.81 | 19.07 | 19.57
ANM+DPSS+SMI | 470.23 | 506.52 | 479.42 | 941.30

To show that the proposed “IVDST+DPSS+SMI” method is suitable for problems
with larger size, we repeat the above experiment with $M=300$,
$(f_{1}^{o},f_{2}^{o},f_{3}^{o})=(0.2,0.24,0.3)$, and $L=2$. Other parameters
are set the same as in the above experiment. Since the SDP-based methods run
slowly in this case, we compare only the “IVDST+DPSS+SMI” method with the SMI
method. We present the time-varying frequency offsets, the dual polynomials,
and the radiation pattern in Figures 5-7. It can be seen that proposed
“IVDST+DPSS+SMI” method still works very well in this high-dimensional case
and significantly outperforms SMI.

(a) Static frequency offset

(b) Linear frequency offset

(c) Zigzag frequency offset

(d) Random frequency offset

Figure 5: Four types of frequency offsets used in the experiments. ($M=300$)

(a) Static frequency offset

(b) Linear frequency offset

(c) Zigzag frequency offset

(d) Random frequency offset

Figure 6: Dual polynomial obtained from the three methods in the four types of
frequency offsets. ($M=300$)

(a) Static frequency offset

(b) Linear frequency offset

(c) Zigzag frequency offset

(d) Random frequency offset

Figure 7: Interference cancellation with DBF in the four types of frequency
offsets. ($M=300$)

Next, we compare the proposed “2D-ANM+DPSS+SMI” method with the SMI method. We
repeat the first experiment with $M=15$333For zigzag frequency offset, we set
$M=30$. and $(f_{1}^{o},f_{2}^{o},f_{3}^{o})=(0.2,0.7,0.7)$. Other parameters
are set same as the first experiment. We present the four types of frequency
offsets in Figure 8 and the four values of $L$ in Table I. We present the dual
polynomials obtained from the “2D-ANM+DPSS+SMI” method in Figure 9 and the
corresponding radiation pattern in Figure 10. It can be seen that the
“2D-ANM+DPSS+SMI” method still significantly outperforms the SMI method.

(a) Static frequency offset

(b) Linear frequency offset

(c) Zigzag frequency offset

(d) Random frequency offset

Figure 8: Four types of frequency offsets used in the experiments. ($M=15$ or
30)

(a) Static frequency offset

(b) Linear frequency offset

(c) Zigzag frequency offset

(d) Random frequency offset

Figure 9: Dual polynomials obtained from the “2D-ANM+DPSS+SMI” method in the
four types of frequency offsets. The red stars denote the true frequencies and
angles in the desired signal and interferers.

(a) Static frequency offset

(b) Linear frequency offset

(c) Zigzag frequency offset

(d) Random frequency offset

Figure 10: Interference cancellation with DBF in the four types of frequency
offsets.

Finally, we repeat the experiments from the previous paragraph (comparing
“2D-ANM+DPSS+SMI” with SMI) for 20 trials each and present the histogram of
the radiation pattern (dB) evaluated at the two interferer directions
($-60^{\circ}$ and $20^{\circ}$) in Figure 11. In each trial, we generate the
four types of frequencies as follows: (a) Static frequency offset: the values
are set as three random integers between 1 and 6 scaled by a factor of 0.01 or
$-0.01$. (b) Linear frequency offset: the three slopes are set as three random
integers between 1 and 6 scaled by a factor of 0.001 or $-0.001$. (c) Zigzag
frequency offset: the three slopes are set as three random integers between 1
and 6 scaled by a factor of 0.001 or $-0.001$. (d) Random frequency offset: an
independent random copy of the offsets shown in Figure 8 (d). The other
settings are same as above. It can be seen from Figure 11 that in most cases
the proposed “2D-ANM+DPSS+SMI” method can achieve a lower radiation pattern at
both directions $-60^{\circ}$ and $20^{\circ}$ with high probability.

(a) Static frequency offset

(b) Linear frequency offset

(c) Zigzag frequency offset

(d) Random frequency offset

Figure 11: The histogram plots of radiation pattern (in dB) evaluated at
$-60^{\circ}$ and $20^{\circ}$ in the four types of frequency offsets with 20
trials.

##  VIII Conclusion

In this paper, we present two novel ANM-based methods in both 1D and 2D
frameworks to address the dynamic problem of interference cancellation with
time-varying frequency offset. By noting that solving the corresponding SDP is
a computational burden in high-dimensional dataset, we also present a novel
fast algorithm to approximately solve the 1D ANM optimization problem.
Finally, we conduct a series of experiments to confirm the benefits of our 1D
and 2D ANM frameworks compared to the conventional SMI, and we illustrate the
computational speedup provided by our fast 1D algorithm. We leave the work of
developing a fast algorithm for the 2D ANM method as our future work.

## Acknowledgement

MW and SL were supported by NSF grant CCF-1704204.

## References

  * [1] Z. Li, Y. Liu, K. G. Shin, J. Liu, and Z. Yan, “Interference steering to manage interference in iot,” IEEE Internet of Things Journal, vol. 6, no. 6, pp. 10458–10471, 2019.
  * [2] L. Chettri and R. Bera, “A comprehensive survey on Internet of Things (IoT) toward 5G wireless systems,” IEEE Internet of Things Journal, vol. 7, no. 1, pp. 16–32, 2019.
  * [3] B. Widrow, P. Mantey, L. Griffiths, and B. Goode, “Adaptive antenna systems,” Proceedings of the IEEE, vol. 55, no. 12, pp. 2143–2159, 1967.
  * [4] H. L. Van Trees, Optimum array processing: Part IV of detection, estimation, and modulation theory.  John Wiley & Sons, 2004.
  * [5] R. J. Mailloux, Phased array antenna handbook.  Artech house, 2017.
  * [6] B. Yang, Z. Yu, J. Lan, R. Zhang, J. Zhou, and W. Hong, “Digital beamforming-based massive MIMO transceiver for 5G millimeter-wave communications,” IEEE Transactions on Microwave Theory and Techniques, vol. 66, no. 7, pp. 3403–3418, 2018.
  * [7] W. Roh, J.-Y. Seol, J. Park, B. Lee, J. Lee, Y. Kim, J. Cho, K. Cheun, and F. Aryanfar, “Millimeter-wave beamforming as an enabling technology for 5G cellular communications: Theoretical feasibility and prototype results,” IEEE Communications Magazine, vol. 52, no. 2, pp. 106–113, 2014.
  * [8] D. Gaydos, P. Nayeri, and R. Haupt, “Experimental demonstration of a software-defined-radio adaptive beamformer,” in 2018 15th European Radar Conference (EuRAD), pp. 561–564, IEEE, 2018.
  * [9] B. Monzingo, R. Haupt, and T. Miller, Introduction to adaptive arrays.  The Institution of Engineering and Technology, 2011.
  * [10] A.-J. Van Der Veen, “Algebraic methods for deterministic blind beamforming,” Proceedings of the IEEE, vol. 86, no. 10, pp. 1987–2008, 1998.
  * [11] Q. Wu and K. M. Wong, “Blind adaptive beamforming for cyclostationary signals,” IEEE Transactions on Signal Processing, vol. 44, no. 11, pp. 2757–2767, 1996.
  * [12] E. Gonen and J. M. Mendel, “Applications of cumulants to array processing. III. Blind beamforming for coherent signals,” IEEE Transactions on Signal Processing, vol. 45, no. 9, pp. 2252–2264, 1997.
  * [13] C. Coviello and L. Sibul, “Blind source separation and beamforming: algebraic technique analysis,” IEEE Transactions on Aerospace and Electronic Systems, vol. 40, no. 1, pp. 221–235, 2004.
  * [14] J. Liu, W. Liu, H. Liu, B. Chen, X.-G. Xia, and F. Dai, “Average SINR calculation of a persymmetric sample matrix inversion beamformer,” IEEE Transactions on Signal Processing, vol. 64, no. 8, pp. 2135–2145, 2015.
  * [15] D. Gaydos, P. Nayeri, and R. Haupt, “Adaptive beamforming in high-interference environments using a software-defined radio array,” in 2019 IEEE International Symposium on Antennas and Propagation and USNC-URSI Radio Science Meeting, pp. 1501–1502, IEEE, 2019.
  * [16] M. W. Ganz, S. L. Wilson, and R. L. Moses, “Convergence of the SMI and the diagonally loaded SMI algorithms with weak interference,” IEEE Transactions on Antennas and Propagation, vol. 38, pp. 394–399, 1990.
  * [17] R. L. Dilsavor and R. L. Moses, “Analysis of modified SMI method for adaptive array weight control,” IEEE Transactions on Signal Processing, vol. 41, no. 2, pp. 721–726, 1993.
  * [18] J. Gu, “Robust beamforming based on variable loading,” Electronics Letters, vol. 41, no. 2, pp. 55–56, 2005.
  * [19] X. Li, D.-W. Wang, X. Ma, and Z. Xiong, “Robust adaptive beamforming using iterative variable loaded sample matrix inverse,” Electronics Letters, vol. 54, no. 9, pp. 546–548, 2018.
  * [20] S. Li, D. Gaydos, P. Nayeri, and M. B. Wakin, “Adaptive interference cancellation using atomic norm minimization,” in 2020 International Applied Computational Electromagnetics Society Symposium (ACES), pp. 1–2, IEEE, 2020.
  * [21] S. Li, D. Gaydos, P. Nayeri, and M. B. Wakin, “Adaptive interference cancellation using atomic norm minimization and denoising,” IEEE Antennas and Wireless Propagation Letters, vol. 19, no. 12, pp. 2349–2353, 2020\.
  * [22] E. J. Candès et al., “Compressive sampling,” in Proceedings of the International Congress of Mathematicians, vol. 3, pp. 1433–1452, Madrid, Spain, 2006.
  * [23] E. J. Candès and M. B. Wakin, “An introduction to compressive sampling,” IEEE Signal Processing Magazine, vol. 25, no. 2, pp. 21–30, 2008.
  * [24] G. Tang, B. N. Bhaskar, P. Shah, and B. Recht, “Compressed sensing off the grid,” IEEE Transactions on Information Theory, vol. 59, no. 11, pp. 7465–7490, 2013.
  * [25] Y. Li and Y. Chi, “Off-the-grid line spectrum denoising and estimation with multiple measurement vectors,” IEEE Transactions on Signal Processing, vol. 64, no. 5, pp. 1257–1269, 2015.
  * [26] Z. Yang and L. Xie, “Exact joint sparse frequency recovery via optimization methods,” IEEE Transactions on Signal Processing, vol. 64, no. 19, pp. 5145–5157, 2016.
  * [27] D. Slepian, “Prolate spheroidal wave functions, Fourier analysis, and uncertainty—V: The discrete case,” Bell System Technical Journal, vol. 57, no. 5, pp. 1371–1430, 1978.
  * [28] M. A. Davenport and M. B. Wakin, “Compressive sensing of analog signals using discrete prolate spheroidal sequences,” Applied and Computational Harmonic Analysis, vol. 33, no. 3, pp. 438–472, 2012.
  * [29] Z. Zhu and M. B. Wakin, “Approximating sampled sinusoids and multiband signals using multiband modulated DPSS dictionaries,” Journal of Fourier Analysis and Applications, vol. 23, no. 6, pp. 1263–1310, 2017.
  * [30] V. Chandrasekaran, B. Recht, P. A. Parrilo, and A. S. Willsky, “The convex geometry of linear inverse problems,” Foundations of Computational mathematics, vol. 12, no. 6, pp. 805–849, 2012.
  * [31] E. J. Candè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.
  * [32] S. Li, D. Yang, G. Tang, and M. B. Wakin, “Atomic norm minimization for modal analysis from random and compressed samples,” IEEE Transactions on Signal Processing, vol. 66, no. 7, pp. 1817–1831, 2018.
  * [33] Y. Xie, S. Li, G. Tang, and M. B. Wakin, “Radar signal demixing via convex optimization,” in 2017 22nd International Conference on Digital Signal Processing (DSP), pp. 1–5, IEEE, 2017.
  * [34] S. Li, M. B. Wakin, and G. Tang, “Atomic norm denoising for complex exponentials with unknown waveform modulations,” IEEE Transactions on Information Theory, vol. 66, no. 6, pp. 3893–3913, 2020.
  * [35] M. B. Wakin, “A study of the temporal bandwidth of video and its implications in compressive sensing,” Colorado School of Mines Technical Report, pp. 1–50, 2012.
  * [36] J. Helland, M. B. Wakin, and G. Tang, “A super-resolution algorithm for extended target localization,” in 2019 IEEE 8th International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP), pp. 386–390, IEEE, 2019.
  * [37] M. Grant, S. Boyd, and Y. Ye, “CVX: Matlab software for disciplined convex programming,” 2008.
  * [38] R. Tütüncü, K. Toh, and M. Todd, “SDPT3—a Matlab software package for semidefinite-quadratic-linear programming, version 3.0,” Web page http://www. math. nus. edu. sg/mattohkc/sdpt3. html, 2001.
  * [39] B. N. Bhaskar, G. Tang, and B. Recht, “Atomic norm denoising with applications to line spectral estimation,” IEEE Transactions on Signal Processing, vol. 61, no. 23, pp. 5987–5999, 2013.
  * [40] D. P. Bertsekas and J. N. Tsitsiklis, Parallel and distributed computation: numerical methods, vol. 23.  Prentice Hall Englewood Cliffs, NJ, 1989.
  * [41] S. Boyd, N. Parikh, and E. Chu, Distributed optimization and statistical learning via the alternating direction method of multipliers.  Now Publishers Inc, 2011.
  * [42] T. L. Hansen and T. L. Jensen, “A fast interior-point method for atomic norm soft thresholding,” Signal Processing, vol. 165, pp. 7–19, 2019.
  * [43] Y. Wang and Z. Tian, “IVDST: A fast algorithm for atomic norm minimization in line spectral estimation,” IEEE Signal Processing Letters, vol. 25, no. 11, pp. 1715–1719, 2018.
  * [44] Y. Liu, Z. Chu, and Y. Yang, “Iterative vandermonde decomposition and shrinkage-thresholding based two-dimensional grid-free compressive beamforming,” The Journal of the Acoustical Society of America, vol. 148, no. 3, pp. EL301–EL306, 2020.


are indexed on the later metric).

  * •

for $f(M_{2},c_{0}^{\prime},c_{1}^{\prime},\infty,s,d)$, we look for the last
element in $\texttt{pfront}[c_{0}^{\prime}-1][d]$ verifying constraint
$c_{1}^{\prime}$. The path minimizing $M_{2}$ being, by design, the last
element.

  * •

to compute $f(M_{0},\infty,c_{1}^{\prime},c_{2}^{\prime},s,d)$, let us first
denote $k$ the smallest integer such that $\texttt{pfront}[k][d]$ contains an
element verifying constraints $c_{1}^{\prime}$ and $c_{2}^{\prime}$. The
resulting image is then any of such elements in $\texttt{pfront}[k][d]$.

As one might notice, computing
$f(M_{j},c_{0}^{\prime},c_{1}^{\prime},c_{2}^{\prime},s,d)$, $j=0,1,2$ cannot
always be achieved in constant time (for $j=0$ and sub-constraints in
particular). Indeed, we favor a simple data structure. A search in an ordered
list of size $\Gamma$ is needed for stricter constraints (and may be performed
$\log(MSD)$ times when optimizing $M_{0}$). To improve the time efficiency of
our solution, each $\texttt{pfront}[i][d]$ may be defined as or converted into
a static array in the implementation.

Finally, for simplicity, we did not show in our pseudo-code the structure and
operations that store and extend the lists of segments. In practice, we store
one representative of the best predecessors and a posteriori retrieve the
lists using a backward induction for each destination.

##  Appendix B SR Graph & Live Conversion Algorithm

The multi-metric SR Graph may be used in different ways. The path computation
algorithm may be run directly on it. While this induces the cost of exploring
a fully-meshed graph, it allows to treat the segment metric as a standard
graph metric (in the form of the number of edges). Another way to use the SR
Graph is to use the information it contains to perform conversions. The
algorithm thus explores the original, sparser topology and converts the
explored paths into segment lists on the fly to consider their number of
segments. Note that the information within the SR Graph may be stored in
different fashions (_e.g.,_ by keeping the results of the APSP computation
separated by source rather than merging it in a single SR Graph). However, the
multi-metric APSP computation itself is mandatory.

The second method requires to design an efficient conversion algorithm, able
to transform a path into a minimal segment list which respects the cost _and_
delay of the original path. This conversion is not trivial, as one must take
into account the (forced) path diversity brought by ECMP. Indeed, using a node
segment implies that the packet may follow any of the ECMP paths, which may
possess heterogeneous delays. On the other hand, relying too much on adjacency
segments will not lead to the minimal segment encoding of the considered
paths. In addition, the number of segments, being an "off the graph" metric,
exhibits peculiar behavior and requires to slightly extend the way path are
checked for dominancy.

###  B.1 Conversion Algorithm

We start by introducing the following notations. We consider a network
$G=(V,E)$. Let $p=(x_{0},\dots,x_{l})$ be a path within $G$, with $x_{i}\in
V$. Let $p_{[x_{i},x_{j}]}=(x_{i},\dots,x_{j})$. We note $d(p)$ the couple of
distances $(d_{1}(p),d_{2}(p))$ of the path $p$.

We denote $\mathit{Dag}$ a shortest path tree (regarding the $M_{2}$ metric,
_i.e.,_ the IGP cost). $\mathit{Dag}_{[u,v]}$ denotes the paths and distances
from $u$ to $v$, within the shortest path tree $\mathit{Dag}$ rooted at $u$.
To encode a given path, the $\mathit{Dag}$ rooted at each node $u\in V$ must
be known. The latter must account for all ECMP path. Alternatively, one may
compute the SR graph of the graph $G$, which contains the minimal amount of
information required to perform any path encoding.

We denote $SR_{G}$ the SR graph of $G$, computed as explained in the main body
of the paper (Section 3.2), but with the adjacency segments removed (the
latter are not necessary to perform path encoding). We denote segments
$S_{i}=(t_{i},s_{i},s_{i+1})$ with $s_{i}\in V$ and
$t_{i}\in\\{\mathit{AdjSeg},\mathit{NodeSeg}\\}$. The latter may either encode
the best ECMP paths from $s_{i}$ to $s_{i+1}$ (in which case,
$t_{i}=\mathit{NodeSeg}$) or a specific link between the two nodes (in which
case, $t_{i}=\mathit{AdjSeg}$). We denote $SR(u,v)$ the best guaranteed
distances when using a node segment from $u$ to reach $v$ (_i.e.,_ the maximum
delay among all path within $\mathit{Dag}_{[u,v]}$). In other words, it is the
weight of the edge $(u,v)$ within $SR_{G}$. We define a list of segments as
$S=(S_{i})_{0\leq i\leq l-1}$.

We now define the _encoding_ of a path $p$.

######  Definition 1 (Encoding)

Let $p=(x_{0},\dots,x_{l})$ be a path. Let $S=(S_{i})_{0\leq i\leq l-1}$ be a
segment list. $S$ encodes the path $p$ if both conditions are verified:

  * •

$\forall i$, $\exists j$ such that $s_{i}=x_{j}$

  * •

$\forall i$, either

    * –

$p_{[s_{i},s_{i+1}]}\in Dag_{[s_{i},s_{i+1}]}\wedge SR(s_{i},s_{i+1})\leq
d(p_{[s_{i},s_{i+1}]})$

    * –

$t_{i}=\mathit{AdjSeg}\wedge(s_{i},s_{i+1})\in p$.

We then define the _delay-cost path encoding problem_.

######  Definition 2 (Delay-cost path encoding problem)

Given a path $p$, the delay-cost path encoding problem consists in finding a
segment list $L$ that encodes $p$ with the minimal number of segments.

To solve the delay-cost path encoding problem, we start by presenting an
algorithm which finds, given a path $p$, the longest prefix of $p$ encodable
by a single segment. The pseudo-code of the algorithm is shown in Alg. 4.

1 if $p_{[x_{0},x_{1}]}\notin\mathit{Dag}_{[x_{0},x_{1}]}$ then

2 return $(\mathit{AdjSeg},x_{0},x_{1})$

3

4

5for $i=2\dots l$ do

6 if $p_{[x_{0},x_{i}]}\notin\mathit{Dag}_{[x_{0},x_{i}]}$ or
$SR(x_{0},x_{i})>d(p_{[x_{0},x_{i}]})$ then

7 return $\mathit{NodeSeg},x_{0},x_{i-1})$

8

9

10return $\mathit{NodeSeg},x_{0},x_{l})$

11

Algorithm 4 1SegLongestPrefix($p$))

The test performed at Line 4 checks whether the first edge of $p$ is within
the shortest path tree of rooted at $x_{0}$. If so, the edge is not part of
the ECMP path from $x_{0}$ to $x_{1}$. Thus, an adjacency segment is required.
Otherwise, the algorithms checks, for each $i=2\dots l$ whether the path
$p_{[x_{0},x_{i}]}$ is encodable in a single segment. This is equivalent to
checking whether the path is within the shortest path tree rooted at $x_{0}$,
and if the latter has, among all the ECMP paths
$\mathit{Dag}_{[x_{0},x_{i}]}$, the highest delay. Once this checks does not
hold anymore, the algorithm returns.

Note that because the SR graph $SR_{G}$ encodes the minimal amount of
information for such conversion, both conditions (Line 4 and 4) can be checked
easily and efficiently by relying on the SR graph $SR_{G}$. Indeed, the
condition at Line 4 can be expressed equivalently as
$SR(x_{0},x_{i-1})+d(x_{i-1},x_{i})\neq SR(x_{0},x_{i})$. Similarly, Line 4
can be expressed as $SR(x_{0},x_{1})\neq d(x_{0},x_{1})$.

We then generalize this algorithm to encode any given path $p$ with a minimal
number of segments, as shown in Alg. 5.

1 $s:=x_{0}$

2 $S:=\emptyset$

3 do

4 $S_{i}=(t,s,v)$ = 1SegLongestPrefix($p_{[s,x_{l}]}$)

5 $L$.push($S_{i}$)

6 $s:=v$

7

8while $s\neq x_{l}$;

9return $L$

Algorithm 5 encode($p$)

Alg. 5 starts by finding the segment $S_{i}$ encoding the longest prefix in
$p$ encodable in a single segment. Suppose that $S_{i}$ encodes the subpath
$p_{[s,v]}$. The algorithm adds the segment to the segment list $S$ and
repeats the procedure, considering the reminder of the path $p_{[v,x_{l}]}$
and $v$ as the new source. This is repeated until $v$ is equal to the last
node in $p$, _i.e.,_ $x_{l}$.

###### Lemma 1

Let $S$ be a segment list that encodes $p$. Then, $|S|\geq|$encode($p$)$|$.

Proof: Let us consider two segment lists $S=((S_{i})_{0\leq i\leq l-1})$ and
$S^{\prime}=((S^{\prime}_{i})_{0\leq i\leq l^{\prime}-1})$, both encoding $p$.
Let $S^{\prime}=encode(p)$ (_i.e.,_ $S^{\prime}$ is the result of Alg. 5).

For the sake of contradiction, let us assume that $l<l^{\prime}$. Then, there
must exist $k$, $0\leq k\leq l-1$, such that

  1. 1.

$s_{k}$ appears before, or at the same place as $s^{\prime}_{k}$ in path $p$

  2. 2.

$s_{k+1}$ appears strictly after $s^{\prime}_{k+1}$ in path $p$.

Let us start by considering that $s_{k}=s^{\prime}_{k}$. Then, 2. contradicts
with the fact that $S^{\prime}_{k}$ is a path encoding the longest prefix of
$p$ (recall that $S_{k}^{\prime}$ is the output of procedure 1SegLongestPrefix
on $s_{k}^{\prime}$).

Otherwise, $s_{k}$ appears strictly before $s^{\prime}_{k}$ in path $p$.
Following the definition of encoding,
$p_{[s_{k},s_{k+1}]}\subset\mathit{Dag}_{[s_{k},s_{k+1}]}$. In other words,
$p_{[s_{k},s_{k+1}]}$ is a shortest path (regarding $M_{2}$). Since
$s^{\prime}_{k}\in p_{\left[s_{k},s_{k+1}\right]}$, then,
$p_{\left[s_{k}^{\prime},s_{k+1}\right]}$ is a shortest path as well, in other
words

$p_{\left[s_{k}^{\prime},s_{k+1}\right]}\subset\mathit{Dag}_{[s^{\prime}_{k},s_{k+1}]}$

and by definition of $SR$,

$d(p_{[s^{\prime}_{k},s_{k+1}]})=d(s^{\prime}_{k},s_{k+1})\leq
SR(s^{\prime}_{k},s_{k+1})$

meaning that $(\mathit{NodeSeg},s^{\prime}_{k},s_{k+1})$ encodes
$p_{[s^{\prime}_{k},s_{k+1}]}$ which contradicts the fact that
$S^{\prime}_{k}$ is the longest encoding.

$\square$

Note that within our code, the conversion algorithm is implemented in a
slightly different fashion for performance sake. During the exploration of the
graph, path are extended edges by edges. To avoid recomputing the encoding of
a path from scratch after each extension, we save for each path states
allowing to efficiently update, for a given path $p$, the number of segments
required if the latter is extended by an edge $e$.

###  B.2 Dominancy checks

Although the algorithms presented in the previous section allow to find the
number of segments to encode a path, it is not sufficient to correctly take
into consideration the number of segments when computing paths. Indeed,
conventional DCLC algorithm tend to only extend non-dominated paths. However,
some precautions must be taken when comparing the number of segments of two
paths. Indeed, because of the peculiar nature of the metric, the latter may
evolve differently. As a consequence, paths that seem dominated must sometimes
be extended nevertheless, as they may end up requiring less segments than
their dominating competitors.

To illustrate this issue, let us consider two paths $p$ and $p^{\prime}$. Let
$S=(S_{i})_{0\leq i\leq l-1}=$ encode$(p)$ and
$S^{\prime}=(S^{\prime}_{i})_{0\leq i\leq l^{\prime}-1}=$
encode($p^{\prime}$). Let
$S^{\prime}_{l-1}=(\mathit{AdjSeg},s^{\prime}_{l-1},s^{\prime}_{l})$ and
$S_{l-1}=(\mathit{NodeSeg},s_{l-1},s_{l})$. Then, intuitively, if $p^{\prime}$
is to be extended by an edge $e=(u,v)$, a new segment is mandatory, as
$S^{\prime}_{l-1}$ may not encode more than a single link. However, a new
segment is not necessarily mandatory to extend $p$, as $S_{l-1}$ may very well
be extended to include $e$, if $p_{[s_{l-1},v]}\in\mathit{Dag}_{[s_{l-1},v]}$
and if $SR(s_{l-1},v)<d(p_{[s_{l},v]})$.

Note that this effect may also occur (although less obviously) even when both
$p$ and $p^{\prime}$ end with a node segment. Indeed, it is possible that
$p_{[s_{l-1},v]}\in\mathit{Dag}_{[s_{l-1},v]}$ while
$p_{[s^{\prime}_{l-1},v]}\notin\mathit{Dag}_{[s^{\prime}_{l-1},v]}$

Consequently, even though a path may seem dominated, this may not be
definitive. Thus, some paths that seem dominated must be nevertheless be
extended. More precisely, a path $p$ _should not_ be considered dominated by a
path $p^{\prime}$, regardless of their cost and delay, if the three following
conditions are met:

  1. 1.

$d_{0}(p)=d_{0}(p^{\prime})$

  2. 2.

$s_{l-1}\neq s^{\prime}_{l-1}$

  3. 3.

$t_{l-1}=\mathit{NodeSeg}$

By modifying the dominancy condition as such, and relying on the conversion
algorithms presented, a standard multi-metric path computation algorithm is
able to correctly take into account the number of segments as another path
metric. Although the overhead induced by our conversion algorithm should be
limited, note that this dominancy condition implies that some paths that are
ultimately dominated may have to be extended nonetheless, increasing the
overall worst-case complexity of the algorithm. The optimal usage of the SR
Graph thus depends on the overhead of the conversion and extended dominancy
check, compared to the cost of exploring a fully meshed graph.

##  Appendix C BEST2COP upper bound with multi-threading

Figure 9: Upper bound of BEST2COP execution time regarding the number of nodes
and threads.

To test BEST2COP upper bound, we force it to explore its full iteration space.
In other words, the algorithm behaves as if $|V|^{2}\times L\times\Gamma$
distances have to be extended at each iteration. The results are shown in Fig.
9 for an increasing number of nodes and threads. BEST2COP does not exceed 84s
when using a single thread, considering $|V|=1000$ and $L=2$, the average
number of parallel links in the transformed graph. This time, reasonable given
the unrealistic nature of the experience, is significantly reduced when
relying on multi-threading. Using 8 threads, BEST2COP execution time decreases
to around 10s, highlighting both the parallelizable nature of our algorithm
and its inherently good performance (past this number of threads, there is no
more speedup as we exceed the number of physical cores of our machine).
Additional experiments conducted on a high-performance grid show that BEST2COP
reaches a speedup factor of 23 when run on 30 cores.