# Nearly optimal resolution estimate for the two-dimensional super-resolution
and a new algorithm for direction of arrival estimation with uniform
rectangular array ††thanks: This work was supported in part by the Swiss
National Science Foundation grant number 200021–200307.

Ping Liu Department of Mathematics, ETH Zürich, Rämistrasse 101, CH-8092
Zürich, Switzerland<EMAIL_ADDRESS>habib.ammari@math.ethz.ch).
Habib Ammari22footnotemark: 2

###### Abstract

In this paper, we develop a new technique to obtain nearly optimal estimate of
the computational resolution limit introduced in [44, 43, 42] for two-
dimensional super-resolution problems. Our main contributions are fivefold:
(i) Our work improves the resolution estimate for number detection and
location recovery in two-dimensional super-resolution problems to nearly
optimal; (ii) As a consequence, we derive a stability result for a sparsity-
promoting algorithm in two-dimensional super-resolution problems (or Direction
of Arrival problems (DOA)). The stability result exhibits the optimal
performance of sparsity promoting in solving such problems; (iii) Our
techniques pave the way for improving the estimate for resolution limits in
higher-dimensional super-resolutions to nearly optimal; (iv) Inspired by these
new techniques, we propose a new coordinate-combination-based model order
detection algorithm for two-dimensional DOA estimation and theoretically
demonstrate its optimal performance, and (v) we also propose a new coordinate-
combination-based MUSIC algorithm for super-resolving sources in two-
dimensional DOA estimation. It has excellent performance and enjoys many
advantages compared to the conventional DOA algorithms. The coordinate-
combination idea seems to be a promising way for multi-dimensional DOA
estimation.

Mathematics Subject Classification: 94A08,94A12, 42A05, 65J22, 65F99,65K99

Keywords: two-dimensional super-resolution, direction of arrival algorithms,
resolution estimates, stability results, sparsity-promoting algorithm, model
order detection, MUSIC algorithm

##  1 Introduction

It is well-known that the physical nature of wave propagation and diffraction
imposes a fundamental barrier in the resolution of imaging systems, which is
termed diffraction limit or resolution limit. Since the famous works of Abbe
[1] and Rayleigh [55] for quantifying the resolution limit, it is widely used
in practice to date that the resolution limit is near half of the wavelength
(see, for instance, [5, 6]). Although this kind of resolution limit was widely
used, it is lack of mathematical foundations and not that applicable to modern
imaging modalities [54, 15]. From the mathematical perspective, the resolution
limit could only be set when taking into account the noise [23, 21, 12] and
surpassing these classical resolution limits is very promising for imaging
modalities with high signal-to-noise ratio (SNR). This understanding motivates
new works on deriving more rigorous resolution limits [30, 31, 48, 47]. At the
beginning of this century, the dependence of two-point resolution on the noise
level has been thoroughly investigated from the perspective of statistical
inference [61, 62, 63], but the resolution estimates for resolving multiple
sources only achieve breakthroughs in recent years due to its nonlinearity.

To understand the resolution in resolving multiple sources, in the earlier
works [44, 43, 42] we have defined “computational resolution limits” for
number detection and location recovery in the one- and multi-dimensional
super-resolution problems and characterized them by the signal-to-noise ratio,
cutoff frequency, and number of sources. In [43], we derived sharp estimates
for the computational resolution limits in one dimensional super-resolution
problems. We extended the estimations to multi-dimensional cases in [42], but
the new estimation is not that sharp due to the techniques of projection used
there. Specifically, the upper bound for the resolution increases rapidly as
the source number $n$ and space dimensionality $k$ increases. To address this
issue, this paper aims to derive better and nearly optimal estimates for the
computational resolution limits in two-dimensional super-resolution problems
and provide a better way to tackle general multi-dimensional cases. The main
contribution of our work are fivefold: (i) Our work improves the resolution
estimate in [42] for number detection and location recovery in two-dimensional
super-resolution problems to nearly optimal; (ii) As a consequence, we derive
a stability result for a sparsity-promoting algorithm in two-dimensional
super-resolution problems (or Direction of Arrival problems (DOA)). Although
it is well-known that the total variation optimization [11] and many other
convex optimization based algorithms [68] have a resolution limit near the
Rayleigh limit [67, 16, 22], our stability result exhibits the optimal super-
resolution ability of $l_{0}$-minimization in solving such problems; (iii) Our
techniques reduce the resolution limit problem to a geometric problem, which
paves the way for improving the estimate for resolution limits in higher
dimensions to nearly optimal; (iv) Inspired by the techniques used in the
proofs, we propose a new coordinate-combination-based model order detection
algorithm for two-dimensional DOA problems and demonstrate its optimal
performance both theoretically and numerically, and (v) we also propose a new
coordinate-combination-based MUSIC (states for MUltiple SIgnal Classification)
algorithm for super-resolving sources in two-dimensional DOA estimation. Our
original algorithm enjoys certain advantages compared to the conventional DOA
algorithms. We also exhibit numerically the phase transition phenomenon of the
algorithm, which demonstrates its excellent resolving capacity. The
coordinate-combination idea seems to be a promising direction for multi-
dimensional DOA estimations.

###  1.1 Existing works on the resolution limit problem

The first theory for quantifying the resolution limit was derived by Ernst
Abbe [1, 71]. Since then, there have been various proposals for the resolution
limit [55, 64, 59, 33], among which the famous and widely used ones are the
Rayleigh limit [55] and the full width at half maximum (FWHM) [20]. However,
these classical resolution limits neglect the effect of noise and hence are
not mathematically rigorous [23, 21, 12]. From a mathematical perspective,
there is no resolution limit when one has perfect access of the exact
intensity profile of the diffraction images. Therefore, the resolution limit
can only be rigorously set when taking into account the measurement noise or
aberration to preclude perfect access to the diffraction images. Based on this
understanding, many works were devoted to characterize the dependence of the
two-point resolution on the signal-to-noise ratio from the perspective of
statistical inference [30, 31, 48, 47, 61, 62, 63]. These classical and semi-
classical limits of two-point resolution have been well-studied and we refer
the reader to [41, 12, 17, 21] for more detailed introductions.

For the resolution limit of superresolving multiple point sources, the problem
becomes much more difficult due to the high degree of nonlinearity. To our
knowledge, the first breakthrough was achieved by Donoho in 1992 [24]. He
considered a grid setting where a discrete measure is supported on a lattice
(spacing by $\Delta$) and regularized by a so-called "Rayleigh index" $b$. The
problem is to reconstruct the amplitudes of the grid points from their noisy
Fourier data in $[-\Omega,\Omega]$ with $\Omega$ being the band limit. He
demonstrated that the minimax error for the amplitude reconstruction is
bounded from below and above by $SRF^{2b-1}\sigma$ and $SRF^{2b+1}\sigma$
respectively with $\sigma$ being the noise level and the super-resolution
factor $SRF=1/({\Omega\Delta})$. His results emphasize the importance of
sparsity and signal-to-noise in super-resolution. But the estimate has not
been improved until recent years. In recent years, due to the enormous
development of super-resolution modalities in biological imaging [29, 70, 32,
10, 57] and the popularity of researches of super-resolution algorithms in
applied mathematics [11, 7, 25, 53, 69, 68, 51, 50, 22, 40, 38], the inherent
superresolving capacity of the imaging problem is drawing increasing interest
and has been well-studied for the one-dimensional case. In [19], the authors
considered resolving $n$-sparse point sources supported on a grid and improved
the results of Donoho. They showed that the minimax error in the amplitude
recovery scales as $SRF^{2n-1}\sigma$ in the presence of noise with intensity
$\sigma$. The case of multi-clustered point sources was considered in [37, 8]
and similar minimax error estimations were derived. In [4, 9], the authors
considered the minimax error for recovering off-the-grid point sources. Based
on an analysis of the "prony-type system", they derived bounds for both
amplitude and location reconstructions of the point sources. More precisely,
they showed that for $\sigma\lessapprox(SRF)^{-2p+1}$, where $p$ is the number
of point sources in a cluster, the minimax error for the amplitude and the
location recoveries scale respectively as $(SRF)^{2p-1}\sigma$ and
$(SRF)^{2p-2}{\sigma}/{\Omega}$, while for the isolated non-clustered source,
the corresponding minimax error for the amplitude and the location recoveries
scale respectively as $\sigma$ and ${\sigma}/{\Omega}$. We also refer the
reader to [49, 12] for understanding the resolution limit from the perceptive
of sample complexity and to [67, 16] for the resolving limit of some
algorithms.

On the other hand, in order to characterize the exact resolution rather than
the minimax error in recovering multiple point sources, in the earlier works
[44, 42, 43] we have defined "computational resolution limits" which
characterize the minimum required distance between point sources so that their
number and locations can be stably resolved under certain noise level. By
developing a nonlinear approximation theory in a so-called Vandermonde space,
we have derived sharp bounds for computational resolution limits in the one-
dimensional super-resolution problem. In particular, we have showed in [43]
that the computational resolution limits for the number and location
recoveries should be respectively
$\frac{C_{\mathrm{num}}}{\Omega}(\frac{\sigma}{m_{\min}})^{\frac{1}{2n-2}}$
and
$\frac{C_{\mathrm{supp}}}{\Omega}(\frac{\sigma}{m_{\min}})^{\frac{1}{2n-1}}$,
where $C_{\mathrm{num}}$ and $C_{\mathrm{supp}}$ are constants and $m_{\min}$
is the minimum strength of the point sources. We have extended these estimates
to multi-dimensional cases in [42] but the results are not that optimal due to
the projection techniques used there. In this paper, we improve the estimates
for the two-dimensional super-resolution problem by a new technique. The
improvements shall be discussed in detail in Section 2. Also, our new
technique paves the way for improving the results in higher-dimensional super-
resolution problems.

###  1.2 Direction of Arrival estimation

Our work also inspires new ideas for the two-dimensional direction of arrival
estimation. Direction of arrival (DOA) estimation refers to the process of
retrieving the direction information of several electromagnetic waves/sources
from the received data of a number of antenna elements in a specific array. It
is an important problem in array signal processing and finds wide applications
in radar, sonar, wireless communications, etc; see, for instance, [6].

In one-dimensional DOA estimation, if the antenna elements are uniformly
spaced in a line, the well-known MUSIC, ESPRIT algorithms, and other subspace
methods can resolve the direction of each incident signal/source with high
resolution. But for the two-dimensional DOA estimation with regular
rectangular array (URA) where both azimuth and elevation angles should be
determined, these subspace methods cannot be simply extended to the two-
dimensional case to directly determine the azimuth and elevation angle of each
source. A major idea to solve the two-dimensional DOA problem is to decompose
it into two independent one-dimensional DOA estimations in which the subspaces
methods can be leveraged to efficiently restore the direction components of
sources corresponding to $x$-axis and $y$-axis. We call the methods with this
decoupling idea as one-dimensional-based algorithms throughout the paper for
convenience of discussion. It is worth emphasizing that other ways for
directly obtaining the azimuth and elevation angles of each source were also
considered [78, 76, 39], but the signal processing in a higher dimensional
space damped their computational efficiency.

Although the one-dimensional-based algorithms are usually much more
computationally efficient, they still suffer from some issues: (i) the loss of
distance separation for $x$-axis or $y$-axis components; (ii) pair matching of
the estimated elevation and azimuth angles. For the first issue, the $x$-axis
(or $y$-axis) components of two sources may be closely spaced even though the
two sources are far away in the two-dimensional space. This causes very
unstable reconstruction of the one-dimensional components and the sources.
Most of the researches usually ignored these issues and some papers proposed
different ways to enhance the reconstruction but the proposed methods are
complicated [73, 72]. For example, in [73], the authors utilized Taylor
expansion, subspace projection, and a tree structure to enhance the
reconstruction when the recovered one-dimensional components are unstable. The
second issue is that the pair matching of the estimated elevation and azimuth
angles is very time consuming when dealing with multiple components of
sources. It usually requires a complex process or two-dimensional search [66,
78, 18, 46, 35].

In this paper, we propose a new efficient one-dimensional-based algorithm for
the two-dimensional DOA estimation which solves the above two issues in a
simple way. First, our algorithm employs a new idea named coordinate-
combination to avoid severe loss of distance separation between sources in
certain region; see Section 5.4 for the detailed discussion. On the other
hand, unlike conventional one-dimensional-based algorithms, the pair matching
problem of our algorithm is a simple balanced assignment problem [52] which
can be solved efficiently by many algorithms such as the Hungarian algorithm.

###  1.3 Organization of the paper

The rest of the paper is organized in the following way. In Section 2, we
present the main results on computational resolution limits for the number
detection and the location recovery in the two-dimensional super-resolution
problem. We also provide a stability result for a sparsity promoting
algorithm. In Section 3, we prove the main results in Section 2. Inspired by
the techniques in the proofs, in Section 4 and Section 5 we introduce
respectively the coordinate-combination-based number detection and source
recovery algorithms in two-dimensional DOA estimations. We also conduct
numerical experiments to demonstrate their super-resolution capability.
Section 6 presents a nonlinear approximation theory in Vandermonde space which
is also a main part in proving our main results. Section 7 is devoted to some
conclusions and future works. In the appendix, we prove a technical lemma.

##  2 Main results

###  2.1 Model setting

We consider the following model of a linear combination of point sources in a
two-dimensional space:

$\mu=\sum_{j=1}^{n}a_{j}\delta_{\boldsymbol{\mathbf{y}}_{j}},$

where $\delta$ denotes Dirac’s $\delta$-distribution in $\mathbb{R}^{2}$,
$\boldsymbol{\mathbf{y}}_{j}\in\mathbb{R}^{2},1\leq j\leq n$, which are the
supports of the measure, represent the locations of the point sources and
$a_{j}\in\mathbb{C},1\leq j\leq n$, their amplitudes. We remark that,
throughout the paper, we will use bold symbols for vectors and matrices, and
ordinary ones for scalar values. We call that the measure $\mu$ is $n$-sparse
if all $a_{j}$’s are nonzero. We denote by

$m_{\min}=\min_{j=1,\cdots,n}|a_{j}|,\quad D_{\min}=\min_{p\neq
j}\Big{|}\Big{|}{\boldsymbol{\mathbf{y}}_{p}-\boldsymbol{\mathbf{y}}_{j}}\Big{|}\Big{|}_{1}.$
(2.1)

We assume that the available measurement is the noisy Fourier data of $\mu$ in
a bounded domain, that is,

$\mathbf{Y}(\boldsymbol{\mathbf{\omega}})=\mathcal{F}\mu(\boldsymbol{\mathbf{\omega}})+\mathbf{W}(\boldsymbol{\mathbf{\omega}})=\sum_{j=1}^{n}a_{j}e^{i\boldsymbol{\mathbf{y}}_{j}^{\top}\boldsymbol{\mathbf{\omega}}}+\mathbf{W}(\boldsymbol{\mathbf{\omega}}),\
\boldsymbol{\mathbf{\omega}}\in[0,\Omega]^{2},$ (2.2)

where $\mathcal{F}\mu$ denotes the Fourier transform of $\mu$, $\Omega$ is the
cut-off frequency, and $\mathbf{W}$ is the noise. We assume that

$||\mathbf{W}(\boldsymbol{\mathbf{\omega}})||_{\infty}<\sigma,$

where $\sigma$ is the noise level. We are interested in the resolution limit
for a cluster of tightly spaced point sources. To be more specific, we denote
by

$B_{\delta,\infty}(\boldsymbol{\mathbf{x}}):=\Big{\\{}\mathbf{y}\ \Big{|}\
\mathbf{y}\in\mathbb{R}^{2},\
||\boldsymbol{\mathbf{y}}-\boldsymbol{\mathbf{x}}||_{\infty}<\delta\Big{\\}},$

and assume that $\boldsymbol{\mathbf{y}}_{j}\in
B_{\frac{(n-1)\pi}{6\Omega},\infty}(\boldsymbol{\mathbf{0}}),j=1,\cdots,n$, or
equivalently
$||\boldsymbol{\mathbf{y}}_{j}||_{\infty}<\frac{(n-1)\pi}{6\Omega}$.

The inverse problem we are interested in is to recover the discrete measure
$\mu$ from the above noisy measurement $\mathbf{Y}$.

###  2.2 Computational Resolution Limit for number detection in the two-
dimensional super-resolution problem

In this section, we estimate the super-resolving capacity of the source number
detection in two-dimensional super-resolution problems. To be specific, we
will define and characterize a computational resolution limit for the
corresponding number detection problems. Our main results are built upon
delicate analysis of the $\sigma$-admissible measure defined below.

######  Definition 2.1.

Given a measurement $\mathbf{Y}$, we say that
$\hat{\mu}=\sum_{j=1}^{m}\hat{a}_{j}\delta_{\mathbf{\hat{y}}_{j}},\
\mathbf{\hat{y}}_{j}\in\mathbb{R}^{2}$ is a $\sigma$-admissible discrete
measure of $\mathbf{Y}$ if

$||\mathcal{F}\hat{\mu}(\boldsymbol{\mathbf{\omega}})-\boldsymbol{\mathbf{Y}}(\boldsymbol{\mathbf{\omega}})||_{\infty}<\sigma,\
\text{for all}\ \boldsymbol{\mathbf{\omega}}\in[0,\Omega]^{2}.$

Note that the set of $\sigma$-admissible measures of $\mathbf{Y}$
characterizes all possible solutions to the inverse problem with the given
measurement $\mathbf{Y}$. If all $\sigma$-admissible measures have at least
$n$ supports, then detecting the correct source number is possible, for
example by targeting at the sparsest admissible measures. However, if there
exists one $\sigma$-admissible measure with less than $n$ supports, detecting
the source number $n$ is impossible without additional prior information. This
leads to the following new definition of resolution limit, named computational
resolution limit.

######  Definition 2.2.

The computational resolution limit to the number detection problem in two
dimensions is defined as the smallest nonnegative number $\mathcal{D}_{2,num}$
such that for all $n$-sparse measures
$\sum_{j=1}^{n}a_{j}\delta_{\mathbf{y}_{j}},\boldsymbol{\mathbf{y}}_{j}\in
B_{\frac{(n-1)\pi}{6\Omega},\infty}(\boldsymbol{\mathbf{0}})$ and the
associated measurement $\boldsymbol{\mathbf{Y}}$ in (2.2), if

$\min_{p\neq j}||\mathbf{y}_{j}-\mathbf{y}_{p}||_{1}\geq\mathcal{D}_{2,num},$

then there does not exist any $\sigma$-admissible measure with less than $n$
supports for $\mathbf{Y}$.

The above resolution limit is termed “computational resolution limit” to
distinguish it from the classic Rayleigh limit. Compared to the Rayleigh
limit, the definition of the computational resolution limit is more rigorous
from the mathematical perspective. It is related to the noise, by which it is
more applicable for modern imaging techniques. In [44, 43, 42], the authors
defined similar computational resolution limits and present rigorous
estimations for them. Here by the following theorem, we derive a nearly
optimal estimate to the $\mathcal{D}_{2,num}$, which substantially improves
the estimate in [42] for the two-dimensional case.

######  Theorem 2.1.

Let the measurement $\mathbf{Y}$ in (2.2) be generated by a $n$-sparse measure
$\mu=\sum_{j=1}^{n}a_{j}\delta_{\mathbf{y}_{j}},\boldsymbol{\mathbf{y}}_{j}\in
B_{\frac{(n-1)\pi}{6\Omega},\infty}(\boldsymbol{\mathbf{0}})$. Let $n\geq 2$
and assume that the following separation condition is satisfied

$\min_{p\neq j,1\leq p,j\leq
n}\Big{|}\Big{|}\mathbf{y}_{p}-\mathbf{y}_{j}\Big{|}\Big{|}_{1}\geq\frac{16.6\pi(n-1)}{\Omega}\Big{(}\frac{\sigma}{m_{\min}}\Big{)}^{\frac{1}{2n-2}}.$
(2.3)

Then there does not exist any $\sigma$-admissible measures of $\mathbf{Y}$
with less than $n$ supports.

Theorem 2.1 reveals that when $\min_{p\neq j,1\leq p,j\leq
n}\Big{|}\Big{|}\mathbf{y}_{p}-\mathbf{y}_{j}\Big{|}\Big{|}_{1}\geq\frac{16.6\pi(n-1)}{\Omega}\Big{(}\frac{\sigma}{m_{\min}}\Big{)}^{\frac{1}{2n-2}}$,
recovering exactly the source number $n$ is possible. Compared with the
Rayleigh limit $\frac{c_{2}\pi}{\Omega}$, where $c_{2}$ is a constant, Theorem
2.1 also indicates that resolving the source number in the sub-Rayleigh regime
is theoretically possible if the SNR is sufficiently large.

Moreover, the estimate in Theorem 2.1 substantially improves the result in
[42], where the upper bound estimation for the two-dimensional computational
resolution limit is

$\frac{Cn(n-1)}{\Omega}\Big{(}\frac{\sigma}{m_{\min}}\Big{)}^{\frac{1}{2n-2}}$

with $C$ being an explicit constant. By the techniques of this paper, we also
pave the way for estimating the resolution limit for higher dimensions. It is
indicated that we can demonstrate that the corresponding resolution limit in
the $k$-dimensional super-resolution problem can be bounded above by

$\frac{C_{num}(k)(n-1)}{\Omega}\Big{(}\frac{\sigma}{m_{\min}}\Big{)}^{\frac{1}{2n-2}},$

where $C_{num}(k)$ is a constant determined by the space dimensionality. This
substantially improves the result in [42] that the computational resolution
limit is estimated to be bounded above by

$\frac{4.4\pi e\
(\pi/2)^{k-1}(n(n-1)/\pi)^{\xi(k-1)}}{\Omega}\Big{(}\frac{\sigma}{m_{\min}}\Big{)}^{\frac{1}{2n-2}},$

where $\xi(k)=\sum_{j=1}^{k}\frac{1}{j},\ k\geq 1$. By these new estimates, we
get rid of the exponential dependence of the index $n$ on the dimensionality
$k$.

On the other hand, it is already known from [42] that the computational
resolution limit for the number detection in the $k$-dimensional super-
resolution problem is bounded below by
$\frac{C_{1}}{\Omega}\Big{(}\frac{\sigma}{m_{\min}}\Big{)}^{\frac{1}{2n-2}}$
for some constant $C_{1}$. Thus the $\mathcal{D}_{2,num}$ is bounded by

$\frac{C_{1}}{\Omega}\Big{(}\frac{\sigma}{m_{\min}}\Big{)}^{\frac{1}{2n-2}}\leq\mathcal{D}_{2,num}\leq\frac{C_{2}n}{\Omega}\Big{(}\frac{\sigma}{m_{\min}}\Big{)}^{\frac{1}{2n-2}}.$
(2.4)

This estimate is nearly optimal.

The above estimates further indicate a phase transition phenomenon in the two-
dimensional number detection problem. Specifically, by (2.4) we expect the
presence of a line of slope $2n-2$ in the parameter space
$\log(SRF)-\log(SNR)$ above which the source number can be correctly detected
in each realization. This phenomenon is confirmed exactly by the number
detection algorithm (Algorithm 2) later in Section 4.5 and illustrated in
Figure 4.2.

###  2.3 Computational Resolution Limit for location recovery in the two-
dimensional super-resolution problem

We next present our results on the resolution limit for the location recovery
problem in two-dimensions. We first introduce the following concept of
$\delta$-neighborhood of discrete measures.
Define

$B_{\delta,1}(\boldsymbol{\mathbf{x}}):=\Big{\\{}\mathbf{y}\ \Big{|}\
\mathbf{y}\in\mathbb{R}^{2},\
||\boldsymbol{\mathbf{y}}-\boldsymbol{\mathbf{x}}||_{1}<\delta\Big{\\}}.$

######  Definition 2.3.

Let $\mu=\sum_{j=1}^{n}a_{j}\delta_{\boldsymbol{\mathbf{y}}_{j}}$ be a
$n$-sparse discrete measure in $\mathbb{R}^{2}$ and let $\delta>0$ be such
that the $n$ balls $B_{\delta,1}(\boldsymbol{\mathbf{y}}_{j}),1\leq j\leq n$
are pairwise disjoint. We say that
$\hat{\mu}=\sum_{j=1}^{n}\hat{a}_{j}\delta_{\mathbf{\hat{y}}_{j}}$ is within
$\delta$-neighborhood of $\mu$ if each $\mathbf{\hat{y}}_{j}$ is contained in
one and only one of the $n$ balls
$B_{\delta,1}(\boldsymbol{\mathbf{y}}_{j}),1\leq j\leq n$.

According to the above definition, a measure $\hat{\mu}$ in a
$\delta$-neighborhood of $\mu$ preserves the inner structure of the collection
of point sources. For a stable location (or support of measure) recovery
algorithm, the output should be a measure in some $\delta$-neighborhood of the
underlying sources. Moreover, $\delta$ should tend to zero as the noise level
$\sigma$ tends to zero. We now introduce the computational resolution limit
for the support recovery problem. For ease of exposition, we only consider
measures supported in
$B_{\frac{(2n-1)\pi}{12\Omega},\infty}(\boldsymbol{\mathbf{0}})$, where $n$ is
the source number.

######  Definition 2.4.

The computational resolution limit in the two-dimensional location recovery
problem is defined as the smallest non-negative number $\mathcal{D}_{2,supp}$
so that for any $n$-sparse measure
$\mu=\sum_{j=1}^{n}a_{j}\delta_{\mathbf{y}_{j}},\boldsymbol{\mathbf{y}}_{j}\in
B_{\frac{(2n-1)\pi}{12\Omega},\infty}^{k}(\boldsymbol{\mathbf{0}})$ and the
associated measurement $\boldsymbol{\mathbf{Y}}$ in (2.2), if

$\min_{p\neq j,1\leq p,j\leq
n}\Big{|}\Big{|}{\mathbf{y}_{p}-\mathbf{y}_{j}}\Big{|}\Big{|}_{1}\geq\mathcal{D}_{2,supp},$

then there exists $\delta>0$ such that any $\sigma$-admissible measure of
$\mathbf{Y}$ with $n$ supports in $B_{\frac{(2n-1)\pi}{12\Omega}}(\mathbf{0})$
is within $\delta$-neighbourhood of $\mu$.

We have the following estimate for the upper bound of $\mathcal{D}_{2,supp}$.

######  Theorem 2.2.

Let $n\geq 2$. Let the measurement $\boldsymbol{\mathbf{Y}}$ in (2.2) be
generated by a $n$-sparse measure
$\mu=\sum_{j=1}^{n}a_{j}\delta_{\boldsymbol{\mathbf{y}}_{j}},\boldsymbol{\mathbf{y}}_{j}\in
B_{\frac{(2n-1)\pi}{12\Omega},\infty}(\boldsymbol{\mathbf{0}})$ in the two-
dimensional space. Assume that

$D_{\min}:=\min_{p\neq
j}\Big{|}\Big{|}\mathbf{y}_{p}-\mathbf{y}_{j}\Big{|}\Big{|}_{1}\geq\frac{15.3\pi(n-\frac{1}{2})}{\Omega}\Big{(}\frac{\sigma}{m_{\min}}\Big{)}^{\frac{1}{2n-1}}.$
(2.5)

If $\hat{\mu}=\sum_{j=1}^{n}\hat{a}_{j}\delta_{\mathbf{\hat{y}}_{j}}$
supported on $B_{\frac{(2n-1)\pi}{12\Omega},\infty}(\boldsymbol{\mathbf{0}})$
is a $\sigma$-admissible measure of $\boldsymbol{\mathbf{Y}}$, then
$\hat{\mu}$ is in a $\frac{D_{\min}}{2}$-neighborhood of $\mu$. Moreover,
after reordering the $\mathbf{\hat{y}}_{j}$’s, we have

$\Big{|}\Big{|}\mathbf{\hat{y}}_{j}-\mathbf{y}_{j}\Big{|}\Big{|}_{1}\leq\frac{C(n)}{\Omega}SRF^{2n-2}\frac{\sigma}{m_{\min}},\quad
1\leq j\leq n,$ (2.6)

where $SRF:=\frac{\pi}{D_{\min}\Omega}$ is the super-resolution factor and

$C(n)=\frac{(1+\sqrt{3})^{2n-1}2^{5n-1}(2n-1)^{2n-1}\pi}{3^{2n-0.5}}.$

Theorem 2.1 demonstrates that when $\min_{p\neq j,1\leq p,j\leq
n}\Big{|}\Big{|}\mathbf{y}_{p}-\mathbf{y}_{j}\Big{|}\Big{|}_{1}\geq\frac{15.3\pi(n-1)}{\Omega}\Big{(}\frac{\sigma}{m_{\min}}\Big{)}^{\frac{1}{2n-2}}$,
it is possible to recover stably the source locations. For sufficiently large
SNR, the limit in Theorem 2.1 is less than the Rayleigh limit. This indicates
that super-resolution is possible for two-dimensional imaging problems. Also,
the estimate here is better than the one obtained in [42], which is

$\frac{Cn(n-1)}{\Omega}\Big{(}\frac{\sigma}{m_{\min}}\Big{)}^{\frac{1}{2n-1}},$

with an explicit constant $C$. By the techniques of this paper, we also pave
the way for estimating the resolution limit of location recovery in higher
dimensional super-resolution problems. In fact, the corresponding resolution
limit in the $k$-dimensional super-resolution problem can be bounded above by

$\frac{C_{supp}(k)(n-1)}{\Omega}\Big{(}\frac{\sigma}{m_{\min}}\Big{)}^{\frac{1}{2n-1}},$

where $C_{supp}(k)$ is a constant depending on the space dimensionality. This
substantially improves the estimate,

$\frac{5.88\pi
e4^{k-1}((n+2)(n-1)/2)^{\xi(k-1)}}{\Omega}\Big{(}\frac{\sigma}{m_{\min}}\Big{)}^{\frac{1}{2n-1}},\quad\xi(k)=\sum_{j=1}^{k}\frac{1}{j},\
k\geq 1,$

derived in [42]. Again, one can get rid of the exponential dependence of the
index of $n$ on the dimensionality $k$ by using these new estimates.

It has been already shown in [42] that the computational resolution limit for
the location recovery in the $k$-dimensional super-resolution problem is
bounded below by
$\frac{C_{3}}{\Omega}\Big{(}\frac{\sigma}{m_{\min}}\Big{)}^{\frac{1}{2n-1}}$
for some constant $C_{3}$. Thus the $\mathcal{D}_{2,supp}$ is bounded by

$\frac{C_{3}}{\Omega}\Big{(}\frac{\sigma}{m_{\min}}\Big{)}^{\frac{1}{2n-1}}\leq\mathcal{D}_{2,supp}\leq\frac{C_{4}n}{\Omega}\Big{(}\frac{\sigma}{m_{\min}}\Big{)}^{\frac{1}{2n-1}}.$
(2.7)

This estimate is nearly optimal.

On the other hand, (2.7) indicates a phase transition in the location recovery
problem. From (2.7) we expect that there exists a line of slope $2n-1$ in the
parameter space of $\log SRF-\log SNR$ such that the location recovery is
stable in every point above the line. This is confirmed by Algorithm 4 in
Section 5.4.2 and illustrated in Figure 5.2.

###  2.4 Stability of a sparsity-promoting algorithm

Sparsity-promoting algorithms are popular methods in imaging processing and
many other fields. By the results for resolution limit, we can derive a
stability result for a $l_{0}$-minimization in the two-dimensional super-
resolution problems. We consider the following $l_{0}$-minimization problem:

$\min_{\rho\in\mathcal{O}}||\rho||_{0}\quad\text{subject
to}\quad|\mathcal{F}\rho(\boldsymbol{\mathbf{\omega}})-\boldsymbol{\mathbf{Y}}(\omega)|<\sigma,\quad\boldsymbol{\mathbf{\omega}}\in[0,\Omega]^{2},$
(2.8)

where $||\rho||_{0}$ is the number of Dirac masses representing the discrete
measure $\rho$. As a corollary of Theorems 2.1 and 2.2, we have the following
stability result.

######  Theorem 2.3.

Let $n\geq 2$ and $\sigma\leq m_{\min}$. Let the measurement
$\boldsymbol{\mathbf{Y}}$ in (2.2) be generated by a $n$-sparse measure
$\mu=\sum_{j=1}^{n}a_{j}\delta_{\boldsymbol{\mathbf{y}}_{j}},\boldsymbol{\mathbf{y}}_{j}\in
B_{\frac{(2n-1)\pi}{12\Omega},\infty}(\boldsymbol{\mathbf{0}})$ in the two-
dimensional space. Assume that

$D_{\min}:=\min_{p\neq
j}\Big{|}\Big{|}\mathbf{y}_{p}-\mathbf{y}_{j}\Big{|}\Big{|}_{1}\geq\frac{15.3\pi(n-\frac{1}{2})}{\Omega}\Big{(}\frac{\sigma}{m_{\min}}\Big{)}^{\frac{1}{2n-1}}.$
(2.9)

Let $\mathcal{O}$ in the minimization problem (2.8) be
$B_{\frac{(n-1)\pi}{6\Omega},\infty}(\boldsymbol{\mathbf{0}})$, then the
solution to (2.8) contains exactly $n$ point sources. For any solution
$\hat{\mu}=\sum_{j=1}^{n}\hat{a}_{j}\delta_{\mathbf{\hat{y}}_{j}}$, it is in a
$\frac{D_{\min}}{2}$-neighborhood of $\mu$. Moreover, after reordering the
$\mathbf{\hat{y}}_{j}$’s, we have

$\Big{|}\Big{|}\mathbf{\hat{y}}_{j}-\mathbf{y}_{j}\Big{|}\Big{|}_{1}\leq\frac{C(n)}{\Omega}SRF^{2n-2}\frac{\sigma}{m_{\min}},\quad
1\leq j\leq n,$ (2.10)

where $SRF:=\frac{\pi}{D_{\min}\Omega}$ and

$C(n)=\frac{(1+\sqrt{3})^{2n-1}2^{5n-1}(2n-1)^{2n-1}\pi}{3^{2n-0.5}}.$

Theorem 2.3 reveals that sparsity promoting over admissible solutions could
resolve the source locations to the resolution limit level. It provides an
insight that theoretically sparsity-promoting algorithms would have excellent
performance on the two-dimensional super-resolution problems. Especially,
under the separation condition (2.9), any tractable sparsity-promoting
algorithms (such as total variation minimization algorithms [11]) rendering
the sparsest solution could stably reconstruct all the source locations.

##  3 Proofs of the main results

The idea for proving the main results of the paper is to use some new
techniques to reduce the two-dimensional problem to a one-dimensional case.
The reduction techniques are mainly based on the three crucial observations in
the following subsection. The estimation methods for the one-dimensional
super-resolution problem are based on a nonlinear approximation theory in
Vandermonde space, which we present in Section 6.

###  3.1 Three crucial observations

We here introduce three crucial observations that reduce the two-dimensional
super-resolution problem to its one-dimensional analog, by which we are able
to derive the resolution limit theory of this paper. Our observations also
pave the way for extending the resolution estimates to higher dimensions.
Moreover, they inspire a new direction for the DOA algorithms; see Sections 4
and 5.

Translation invariant:
By the translation invariant we mean that if a measure
$\hat{\mu}=\sum_{j=1}^{q}\hat{a}_{j}\delta_{\boldsymbol{\mathbf{\hat{y}}}_{j}}$
is a $\sigma$-admissible measure for the measurement
$\boldsymbol{\mathbf{Y}}$, then
$\hat{\mu}=\sum_{j=1}^{q}\hat{a}_{j}\delta_{\boldsymbol{\mathbf{\hat{y}}}_{j}+\boldsymbol{\mathbf{v}}}$
is a $\sigma$-admissible measure for measurement
$e^{i\boldsymbol{\mathbf{v}}^{\top}\boldsymbol{\mathbf{\omega}}}\boldsymbol{\mathbf{Y}}(\boldsymbol{\mathbf{\omega}})$
for any vector $\boldsymbol{\mathbf{v}}\in\mathbb{R}^{2}$. More precisely, we
have

$\Big{|}{\sum_{j=1}^{q}\hat{a}_{j}e^{i(\boldsymbol{\mathbf{\hat{y}}}_{j}+\boldsymbol{\mathbf{v}})^{\top}\boldsymbol{\mathbf{\omega}}}-e^{i\boldsymbol{\mathbf{v}}^{\top}\boldsymbol{\mathbf{\omega}}}\boldsymbol{\mathbf{Y}}(\boldsymbol{\mathbf{\omega}})}\Big{|}=\Big{|}{\sum_{j=1}^{q}\hat{a}_{j}e^{i\boldsymbol{\mathbf{\hat{y}}}_{j}^{\top}\boldsymbol{\mathbf{\omega}}}-\boldsymbol{\mathbf{Y}}(\boldsymbol{\mathbf{\omega}})}\Big{|}<\sigma,\quad\boldsymbol{\mathbf{\omega}}\in[0,\Omega]^{2}.$
(3.1)

In addition, if for certain $\delta\geq 0$,

$\Big{|}{\sum_{j=1}^{q}\hat{a}_{j}e^{i\boldsymbol{\mathbf{\hat{y}}}_{j}^{\top}\boldsymbol{\mathbf{\omega}}}-\sum_{j=1}^{n}a_{j}e^{i\boldsymbol{\mathbf{y}}_{j}^{\top}\boldsymbol{\mathbf{\omega}}}}\Big{|}<\delta,\quad\boldsymbol{\mathbf{\omega}}\in[0,\Omega]^{2},$
(3.2)

then for any vector $\boldsymbol{\mathbf{v}}\in\mathbb{R}^{2}$,

$\Big{|}{\sum_{j=1}^{q}\hat{a}_{j}e^{i(\boldsymbol{\mathbf{\hat{y}}}_{j}+\boldsymbol{\mathbf{v}})^{\top}\boldsymbol{\mathbf{\omega}}}-\sum_{j=1}^{n}a_{j}e^{i(\boldsymbol{\mathbf{y}}_{j}+\boldsymbol{\mathbf{v}})^{\top}\boldsymbol{\mathbf{\omega}}}}\Big{|}<\delta,\quad\boldsymbol{\mathbf{\omega}}\in[0,\Omega]^{2}.$

Combination of coordinates:
The second observation is that if we suppose that (3.2) is satisfied, we have
a similar estimate for the summation of combinations of
$e^{i\tau\boldsymbol{\mathbf{\hat{y}}}_{j,1}},e^{i\tau\boldsymbol{\mathbf{\hat{y}}}_{j,2}}$
and
$e^{i\tau\boldsymbol{\mathbf{y}}_{j,1}},e^{i\tau\boldsymbol{\mathbf{y}}_{j,2}}$
for certain $\tau$. Specifically, we have the following lemma.

######  Lemma 3.1.

For any integer $t\geq 0$ and $\tau\leq\frac{\Omega}{t}$, the measurement
constraint (3.2) implies

$\Big{|}{\sum_{j=1}^{q}\hat{a}_{j}(e^{ir_{1}}e^{i\tau\boldsymbol{\mathbf{\hat{y}}}_{j,1}}+e^{ir_{2}}e^{i\tau\boldsymbol{\mathbf{\hat{y}}}_{j,2}})^{t}-\sum_{j=1}^{n}a_{j}(e^{ir_{1}}e^{i\tau\boldsymbol{\mathbf{y}}_{j,1}}+e^{ir_{2}}e^{i\tau\boldsymbol{\mathbf{y}}_{j,2}})^{t}}\Big{|}<2^{t}\delta,\quad
r_{1},r_{2}\in\mathbb{R}.$

###### Proof.

Let
$\hat{d}_{j}=e^{ir_{1}}e^{i\tau\boldsymbol{\mathbf{\hat{y}}}_{j,1}}+e^{ir_{2}}e^{i\tau\boldsymbol{\mathbf{\hat{y}}}_{j,2}}$
and
$d_{j}=e^{ir_{1}}e^{i\tau\boldsymbol{\mathbf{y}}_{j,1}}+e^{ir_{2}}e^{i\tau\boldsymbol{\mathbf{y}}_{j,2}}$.
We have

$\displaystyle\Big{|}{\sum_{j=1}^{q}\hat{a}_{j}\hat{d}_{j}^{t}-\sum_{j=1}^{n}a_{j}d_{j}^{t}}\Big{|}=\Big{|}{\sum_{j=1}^{q}\hat{a}_{j}(e^{ir_{1}}e^{i\tau\boldsymbol{\mathbf{\hat{y}}}_{j,1}}+e^{ir_{2}}e^{i\tau\boldsymbol{\mathbf{\hat{y}}}_{j,2}})^{t}-\sum_{j=1}^{n}a_{j}(e^{ir_{1}}e^{i\tau\boldsymbol{\mathbf{y}}_{j,1}}+e^{ir_{2}}e^{i\tau\boldsymbol{\mathbf{y}}_{j,2}})^{t}}\Big{|}$
$\displaystyle=$ $\displaystyle\Big{|}{\sum_{t_{1}+t_{2}=t,0\leq
t_{1},t_{2}\leq t}{t\choose
t_{1}}\Big{(}\sum_{j=1}^{q}\hat{a}_{j}e^{ir_{1}t_{1}}e^{ir_{2}t_{2}}e^{i\tau\boldsymbol{\mathbf{\hat{y}}}_{j,1}t_{1}}e^{i\tau\boldsymbol{\mathbf{\hat{y}}}_{j,2}t_{2}}-\sum_{j=1}^{n}a_{j}e^{ir_{1}t_{1}}e^{ir_{2}t_{2}}e^{i\tau\boldsymbol{\mathbf{y}}_{j,1}t_{1}}e^{i\tau\boldsymbol{\mathbf{y}}_{j,2}t_{2}}\Big{)}}\Big{|}$
$\displaystyle\leq$ $\displaystyle\sum_{t_{1}+t_{2}=t,0\leq t_{1},t_{2}\leq
t}{t\choose
t_{1}}\Big{|}{\sum_{j=1}^{q}\hat{a}_{j}e^{i\tau\boldsymbol{\mathbf{\hat{y}}}_{j,1}t_{1}}e^{i\tau\boldsymbol{\mathbf{\hat{y}}}_{j,2}t_{2}}-\sum_{j=1}^{n}a_{j}e^{i\tau\boldsymbol{\mathbf{y}}_{j,1}t_{1}}e^{i\tau\boldsymbol{\mathbf{y}}_{j,2}t_{2}}}\Big{|}$
$\displaystyle=$ $\displaystyle\sum_{t_{1}+t_{2}=t,0\leq t_{1},t_{2}\leq
t}{t\choose
t_{1}}\Big{|}{\sum_{j=1}^{q}\hat{a}_{j}e^{i(t_{1}\tau,t_{2}\tau)\boldsymbol{\mathbf{\hat{y}}}_{j}}-\sum_{j=1}^{n}a_{j}e^{i(t_{1}\tau,t_{2}\tau)\boldsymbol{\mathbf{y}}_{j}}}\Big{|}$
$\displaystyle<$ $\displaystyle\sum_{t_{1}+t_{2}=t,0\leq t_{1},t_{2}\leq
t}{t\choose t_{1}}\delta\quad\Big{(}\text{by $\tau\leq\frac{\Omega}{t}$ and
(\ref{equ:measureconstraint2})}\Big{)}$ $\displaystyle=$ $\displaystyle
2^{t}\delta.$

∎

This is the key observation of the paper. It reduces the two-dimensional
super-resolution problem to nearly a one-dimensional super-resolution one.
Since it is about the difference between summation of combinations of
$e^{i\tau\boldsymbol{\mathbf{\hat{y}}}_{j,1}},e^{i\tau\boldsymbol{\mathbf{\hat{y}}}_{j,2}}$
and
$e^{i\tau\boldsymbol{\mathbf{y}}_{j,1}},e^{i\tau\boldsymbol{\mathbf{y}}_{j,2}}$,
we refer to this observation as combination of coordinates and call the
elements
$e^{i\tau\boldsymbol{\mathbf{y}}_{j,1}}+e^{i\tau\boldsymbol{\mathbf{y}}_{j,2}}$
coordinate-combined elements. This coordinate-combination technique will be
used in deriving new algorithms for the DOA problem in Sections 4 and 5.

Compared to the projection techniques in [42, 12] which utilize the
measurement constraint only in several one-dimensional spaces to derive
stability results, our formulation utilizes more measurement constraints and
consequently yields better estimates.

Preservation of the separation distance for the coordinate-combined elements:
The last observation is that, for $\boldsymbol{\mathbf{\theta}}_{j}$’s in
$[0,\frac{2\pi}{3}]^{2}$, the coordinate-combined elements
$e^{i\boldsymbol{\mathbf{\theta}}_{j,1}}+e^{i\boldsymbol{\mathbf{\theta}}_{j,2}}$
still preserve the separation distance between the
$\boldsymbol{\mathbf{\theta}}_{j}$’s. This is revealed by Lemma 3.2. Note that
the projection trick in [42, 12] and many conventional two-dimensional DOA
algorithms do not preserve the separation distance between the original
source. This causes many issues in the reconstruction and resolution
estimation. Lemma 3.2 is the main result of this paper by which we could
overcome the above issues and hence find a new way to solve two-dimensional
DOA problems.

######  Lemma 3.2.

For two different vectors
$\boldsymbol{\mathbf{\theta}}_{j}\in\mathbb{[}0,\frac{2\pi}{3}]^{2},j=1,2$
with
$\frac{\pi}{3}\leq\boldsymbol{\mathbf{\theta}}_{j,2}-\boldsymbol{\mathbf{\theta}}_{j,1}\leq\frac{2}{3}\pi,j=1,2$,
if
$||\boldsymbol{\mathbf{\theta}}_{1}-\boldsymbol{\mathbf{\theta}}_{2}||_{1}\geq\Delta$,
then

$\Big{|}{e^{i\boldsymbol{\mathbf{\theta}}_{1,1}}+e^{i\boldsymbol{\mathbf{\theta}}_{1,2}}-(e^{i\boldsymbol{\mathbf{\theta}}_{2,1}}+e^{i\boldsymbol{\mathbf{\theta}}_{2,2}})}\Big{|}\geq\frac{3}{2\pi}\Delta.$

###### Proof.

Note that
$0\leq\boldsymbol{\mathbf{\theta}}_{j,1}<\boldsymbol{\mathbf{\theta}}_{j,2}\leq\frac{2\pi}{3},j=1,2$.
We prove the lemma by considering the following two cases.
Case 1:
$0\leq\boldsymbol{\mathbf{\theta}}_{1,1}\leq\boldsymbol{\mathbf{\theta}}_{2,1}<\boldsymbol{\mathbf{\theta}}_{2,2}\leq\boldsymbol{\mathbf{\theta}}_{1,2}\leq\frac{2\pi}{3}$.
In this case,

$\displaystyle\Big{|}{e^{i\boldsymbol{\mathbf{\theta}}_{1,1}}+e^{i\boldsymbol{\mathbf{\theta}}_{1,2}}-(e^{i\boldsymbol{\mathbf{\theta}}_{2,1}}+e^{i\boldsymbol{\mathbf{\theta}}_{2,2}})}\Big{|}\geq$
$\displaystyle\Big{|}{e^{i\boldsymbol{\mathbf{\theta}}_{2,1}}+e^{i\boldsymbol{\mathbf{\theta}}_{2,2}}}\Big{|}-\Big{|}{e^{i\boldsymbol{\mathbf{\theta}}_{1,1}}+e^{i\boldsymbol{\mathbf{\theta}}_{1,2}}}\Big{|}$
$\displaystyle=$ $\displaystyle
2\Big{(}\cos(\frac{\phi_{2}}{2})-\cos(\frac{\phi_{1}}{2})\Big{)},$

where
$\phi_{j}=\boldsymbol{\mathbf{\theta}}_{j,2}-\boldsymbol{\mathbf{\theta}}_{j,1},j=1,2$.
By the assumption made in the lemma, we have
$\Delta\leq\phi_{1}-\phi_{2}\leq\frac{\pi}{3}$. Note also that
$\frac{\pi}{6}\leq\frac{\phi_{1}+\phi_{2}}{4}\leq\frac{\pi}{3}$. Thus

$2\Big{(}\cos(\frac{\phi_{2}}{2})-\cos(\frac{\phi_{1}}{2})\Big{)}=4\sin(\frac{\phi_{1}+\phi_{2}}{4})\sin(\frac{\phi_{1}-\phi_{2}}{4})\geq
4\sin(\frac{\pi}{6})\sin(\frac{\Delta}{4})\geq\frac{3\Delta}{2\pi},$

where the last inequality uses
$\sin(\frac{\Delta}{4})\geq\frac{3}{\pi}\frac{\Delta}{4}$ for
$0<\frac{\Delta}{4}\leq\frac{\pi}{12}$.
Case 2:
$0\leq\boldsymbol{\mathbf{\theta}}_{1,1}\leq\boldsymbol{\mathbf{\theta}}_{2,1}\leq\boldsymbol{\mathbf{\theta}}_{1,2}\leq\boldsymbol{\mathbf{\theta}}_{2,2}\leq\frac{2\pi}{3}$.
The idea is to calculate the angle between
$e^{i\boldsymbol{\mathbf{\theta}}_{1,1}}+e^{i\boldsymbol{\mathbf{\theta}}_{1,2}}$
and
$e^{i\boldsymbol{\mathbf{\theta}}_{2,1}}+e^{i\boldsymbol{\mathbf{\theta}}_{2,2}}$.
By simple analysis of the angle relations between
$e^{i\boldsymbol{\mathbf{\theta}}_{1,1}},e^{i\boldsymbol{\mathbf{\theta}}_{1,2}},e^{i\boldsymbol{\mathbf{\theta}}_{2,1}},e^{i\boldsymbol{\mathbf{\theta}}_{2,2}}$,
we obtain that the angle between
$e^{i\boldsymbol{\mathbf{\theta}}_{1,1}}+e^{i\boldsymbol{\mathbf{\theta}}_{1,2}}$
and
$e^{i\boldsymbol{\mathbf{\theta}}_{2,1}}+e^{i\boldsymbol{\mathbf{\theta}}_{2,2}}$
is
$\frac{\boldsymbol{\mathbf{\theta}}_{2,1}-\boldsymbol{\mathbf{\theta}}_{1,1}+\boldsymbol{\mathbf{\theta}}_{2,2}-\boldsymbol{\mathbf{\theta}}_{1,2}}{2}$,
which is larger than $\frac{\Delta}{2}$. Thus

$\Big{|}{e^{i\boldsymbol{\mathbf{\theta}}_{1,1}}+e^{i\boldsymbol{\mathbf{\theta}}_{1,2}}-(e^{i\boldsymbol{\mathbf{\theta}}_{2,1}}+e^{i\boldsymbol{\mathbf{\theta}}_{2,2}})}\Big{|}\geq\max\Big{(}\Big{|}{e^{i\boldsymbol{\mathbf{\theta}}_{1,1}}+e^{i\boldsymbol{\mathbf{\theta}}_{1,2}}}\Big{|},\Big{|}{e^{i\boldsymbol{\mathbf{\theta}}_{2,1}}+e^{i\boldsymbol{\mathbf{\theta}}_{2,2}}}\Big{|}\Big{)}\sin(\frac{\Delta}{2}).$

Since
$\frac{\pi}{3}\leq\boldsymbol{\mathbf{\theta}}_{j,2}-\boldsymbol{\mathbf{\theta}}_{j,1}\leq\frac{2}{3}\pi,j=1,2$,
we have

$\max\Big{(}\Big{|}{e^{i\boldsymbol{\mathbf{\theta}}_{1,1}}+e^{i\boldsymbol{\mathbf{\theta}}_{1,2}}}\Big{|},\Big{|}{e^{i\boldsymbol{\mathbf{\theta}}_{2,1}}+e^{i\boldsymbol{\mathbf{\theta}}_{2,2}}}\Big{|}\Big{)}\geq
1.$

Therefore,

$\Big{|}{e^{i\boldsymbol{\mathbf{\theta}}_{1,1}}+e^{i\boldsymbol{\mathbf{\theta}}_{1,2}}-(e^{i\boldsymbol{\mathbf{\theta}}_{2,1}}+e^{i\boldsymbol{\mathbf{\theta}}_{2,2}})}\Big{|}\geq\sin(\frac{\Delta}{2})\geq\frac{3\Delta}{2\pi},$

where the last inequality uses
$\sin(\frac{\Delta}{2})\geq\frac{3}{\pi}\frac{\Delta}{2}$ for
$0<\frac{\Delta}{2}\leq\frac{\pi}{6}$. ∎

###  3.2 Proof of Theorem 2.1

###### Proof.

The proof of this theorem is by contradiction. Suppose that there exists a
measure
$\hat{\mu}=\sum_{j=1}^{q}\hat{a}_{j}\delta_{\boldsymbol{\mathbf{\hat{y}}}_{j}}$
with $q<n$ which is a $\sigma$-admissible measure of
$\boldsymbol{\mathbf{Y}}$. Then, by the measurement constraint (2.2) and
$|\boldsymbol{\mathbf{W}}(\boldsymbol{\mathbf{\omega}})|<\sigma$, we have

$\Big{|}{\sum_{j=1}^{q}\hat{a}_{j}e^{i\boldsymbol{\mathbf{\hat{y}}}_{j}^{\top}\boldsymbol{\mathbf{\omega}}}-\sum_{j=1}^{n}a_{j}e^{i\boldsymbol{\mathbf{y}}_{j}^{\top}\boldsymbol{\mathbf{\omega}}}}\Big{|}<2\sigma,\quad\boldsymbol{\mathbf{\omega}}\in[0,\Omega]^{2}.$
(3.3)

Since $\boldsymbol{\mathbf{y}}_{j}\in[-\lambda,\lambda]^{2}$ with
$\lambda=\frac{(n-1)\pi}{6\Omega}$, by letting
$\boldsymbol{\mathbf{v}}=(0,6\lambda)^{\top}$, we obtain

$\boldsymbol{\mathbf{y}}_{j}+\boldsymbol{\mathbf{v}}\in[-\lambda,\lambda]\times[5\lambda,7\lambda].$
(3.4)

On the other hand, by (3.3) we also get

$\Big{|}{\sum_{j=1}^{q}\hat{a}_{j}e^{i(\boldsymbol{\mathbf{\hat{y}}}_{j}+\boldsymbol{\mathbf{v}})^{\top}\boldsymbol{\mathbf{\omega}}}-\sum_{j=1}^{n}a_{j}e^{i(\boldsymbol{\mathbf{y}}_{j}+\boldsymbol{\mathbf{v}})^{\top}\boldsymbol{\mathbf{\omega}}}}\Big{|}<2\sigma,\quad\boldsymbol{\mathbf{\omega}}\in[0,\Omega]^{2}.$

Thus with a slight abuse of notation, we still denote those
$\boldsymbol{\mathbf{\hat{y}}}_{j}+\boldsymbol{\mathbf{v}}$ and
$\boldsymbol{\mathbf{y}}_{j}+\boldsymbol{\mathbf{v}}$ by
$\boldsymbol{\mathbf{\hat{y}}}_{j}$, $\boldsymbol{\mathbf{y}}_{j}$
respectively and consider them in the rest of the proof. Note that we have

$\boldsymbol{\mathbf{y}}_{j}\in[-\lambda,\lambda]\times[5\lambda,7\lambda],\quad
j=1,\cdots,n.$

Let $\tau=\frac{\Omega}{2(n-1)}$, together with
$\lambda=\frac{(n-1)\pi}{6\Omega}$, we have
$\tau\boldsymbol{\mathbf{y}}_{j}\in[-\frac{\pi}{12},\frac{\pi}{12}]\times[\frac{5\pi}{12},\frac{7\pi}{12}]$.
This yields

$-\frac{\pi}{12}\leq\tau\boldsymbol{\mathbf{y}}_{j,1}\leq\frac{\pi}{12},\quad\frac{5\pi}{12}\leq\tau\boldsymbol{\mathbf{y}}_{j,2}\leq\frac{7\pi}{12},\quad\frac{\pi}{3}\leq\tau\boldsymbol{\mathbf{y}}_{j,2}-\tau\boldsymbol{\mathbf{y}}_{j,1}\leq\frac{2\pi}{3}.$
(3.5)

On the other hand, let
$\hat{d}_{j}=e^{i\tau\boldsymbol{\mathbf{\hat{y}}}_{j,1}}+e^{i\tau\boldsymbol{\mathbf{\hat{y}}}_{j,2}}$
and
$d_{j}=e^{i\tau\boldsymbol{\mathbf{y}}_{j,1}}+e^{i\tau\boldsymbol{\mathbf{y}}_{j,2}}$.
By Lemma 3.1 and (3.3) we have that

$\Big{|}{\sum_{j=1}^{q}\hat{a}_{j}\hat{d}_{j}^{t}-\sum_{j=1}^{n}a_{j}d_{j}^{t}}\Big{|}<2^{t+1}\sigma,\quad
t=0,1,\cdots,2n-2.$ (3.6)

Let

$\boldsymbol{\mathbf{b}}=\Big{(}\sum_{j=1}^{q}\hat{a}_{j}\hat{d}_{j}^{0}-\sum_{j=1}^{n}a_{j}d_{j}^{0},\quad\sum_{j=1}^{q}\hat{a}_{j}\hat{d}_{j}^{1}-\sum_{j=1}^{n}a_{j}d_{j}^{1},\quad\cdots,\quad\sum_{j=1}^{q}\hat{a}_{j}\hat{d}_{j}^{2n-2}-\sum_{j=1}^{n}a_{j}d_{j}^{2n-2}\Big{)}^{\top}.$

Since (3.5) holds, Lemma 3.2 yields

$d_{\min}:=\min_{p\neq
q}\Big{|}{d_{p}-d_{q}}\Big{|}\geq\frac{3}{2\pi}\min_{p\neq
q}\tau\Big{|}\Big{|}\boldsymbol{\mathbf{y}}_{p}-\boldsymbol{\mathbf{y}}_{q}\Big{|}\Big{|}_{1}>12.4\Big{(}\frac{\sigma}{m_{\min}}\Big{)}^{\frac{1}{2n-2}}>2\sqrt{6(1+\sqrt{3})}\Big{(}\frac{4}{\sqrt{3}}\frac{\sigma}{m_{\min}}\Big{)}^{\frac{1}{2n-2}},$

where the second last inequality is due to the separation condition (2.3). On
the other hand, we have $|\hat{d}_{p}|\leq 2,p=1,\cdots,q$ and
$|d_{j}|\leq\sqrt{3},j=1,\cdots,n$ since (3.5) holds. Thus we can apply
Theorem 6.2 and get

$||\boldsymbol{\mathbf{b}}||_{2}\geq\frac{m_{\min}(d_{\min})^{2n-2}}{(2(1+2)(1+\sqrt{3}))^{(n-1)}}>\frac{4^{n}\sigma}{\sqrt{3}}.$

However, (3.6) implies that
$||\boldsymbol{\mathbf{b}}||_{2}<\frac{4^{n}\sigma}{\sqrt{3}}$, which is a
contradiction. This proves the theorem. ∎

###  3.3 Proof of Theorem 2.2

###### Proof.

Note that $\boldsymbol{\mathbf{y}}_{j},\boldsymbol{\mathbf{\hat{y}}}_{j}$’s
are in $[-\lambda,\lambda]^{2}$ with $\lambda=\frac{(2n-1)\pi}{12\Omega}$ and
$\hat{\mu}=\sum_{j=1}^{n}\hat{a}_{j}\delta_{\boldsymbol{\mathbf{\hat{y}}}_{j}}$
is a $\sigma$-admissible measure of $\boldsymbol{\mathbf{Y}}$. Let
$\tau=\frac{\Omega}{2n-1}$. Similarly to the proof in the above section, we
can construct
$\boldsymbol{\mathbf{x}}_{j}=\boldsymbol{\mathbf{y}}_{j}+\boldsymbol{\mathbf{v}},\boldsymbol{\mathbf{\hat{x}}}_{j}=\boldsymbol{\mathbf{\hat{y}}}_{j}+\boldsymbol{\mathbf{v}}$
so that
$\tau\boldsymbol{\mathbf{\hat{x}}}_{j},\tau\boldsymbol{\mathbf{x}}_{j}\in[-\frac{\pi}{12},\frac{\pi}{12}]\times[\frac{5\pi}{12},\frac{7\pi}{12}]$
and

$\Big{|}{\sum_{j=1}^{n}\hat{a}_{j}e^{i\boldsymbol{\mathbf{\hat{x}}}_{j}^{\top}\boldsymbol{\mathbf{\omega}}}-\sum_{j=1}^{n}a_{j}e^{i\boldsymbol{\mathbf{x}}_{j}^{\top}\boldsymbol{\mathbf{\omega}}}}\Big{|}<2\sigma,\quad\boldsymbol{\mathbf{\omega}}\in[0,\Omega]^{2}.$
(3.7)

Thus we have

$-\frac{\pi}{12}\leq\tau\boldsymbol{\mathbf{x}}_{j,1}\leq\frac{\pi}{12},\quad\frac{5\pi}{12}\leq\tau\boldsymbol{\mathbf{x}}_{j,2}\leq\frac{7\pi}{12},\quad\frac{\pi}{3}\leq\tau\boldsymbol{\mathbf{x}}_{j,2}-\tau\boldsymbol{\mathbf{x}}_{j,1}\leq\frac{2\pi}{3},$
(3.8)
$-\frac{\pi}{12}\leq\tau\boldsymbol{\mathbf{\hat{x}}}_{j,1}\leq\frac{\pi}{12},\quad\frac{5\pi}{12}\leq\tau\boldsymbol{\mathbf{\hat{x}}}_{j,2}\leq\frac{7\pi}{12},\quad\frac{\pi}{3}\leq\tau\boldsymbol{\mathbf{\hat{x}}}_{j,2}-\tau\boldsymbol{\mathbf{\hat{x}}}_{j,1}\leq\frac{2\pi}{3}.$
(3.9)

Moreover, it follows that

$-\frac{\pi}{12}\leq\tau\boldsymbol{\mathbf{x}}_{j,1}\leq\frac{\pi}{12},\quad\frac{-7\pi}{12}\leq\tau\boldsymbol{\mathbf{x}}_{j,2}-\pi\leq\frac{-5\pi}{12},\quad\frac{\pi}{3}\leq\tau\boldsymbol{\mathbf{x}}_{j,1}-(\tau\boldsymbol{\mathbf{x}}_{j,2}-\pi)\leq\frac{2\pi}{3},$
(3.10)
$-\frac{\pi}{12}\leq\tau\boldsymbol{\mathbf{\hat{x}}}_{j,1}\leq\frac{\pi}{12},\quad\frac{-7\pi}{12}\leq\tau\boldsymbol{\mathbf{\hat{x}}}_{j,2}-\pi\leq\frac{-5\pi}{12},\quad\frac{\pi}{3}\leq\tau\boldsymbol{\mathbf{\hat{x}}}_{j,1}-(\tau\boldsymbol{\mathbf{\hat{x}}}_{j,2}-\pi)\leq\frac{2\pi}{3}.$
(3.11)

Let
$\hat{d}_{j}=e^{i\tau\boldsymbol{\mathbf{\hat{x}}}_{j,1}}+e^{i\tau\boldsymbol{\mathbf{\hat{x}}}_{j,2}},d_{j}=e^{i\tau\boldsymbol{\mathbf{x}}_{j,1}}+e^{i\tau\boldsymbol{\mathbf{x}}_{j,2}}$
and
$\hat{g}_{j}=e^{i\tau\boldsymbol{\mathbf{\hat{x}}}_{j,1}}+e^{i(\tau\boldsymbol{\mathbf{\hat{x}}}_{j,2}-\pi)},g_{j}=e^{i\tau\boldsymbol{\mathbf{x}}_{j,1}}+e^{i(\tau\boldsymbol{\mathbf{x}}_{j,2}-\pi)}$.
By (3.7) and Lemma 3.1, we arrive at

$\displaystyle\Big{|}{\sum_{j=1}^{n}\hat{a}_{j}\hat{d}_{j}^{t}-\sum_{j=1}^{n}a_{j}d_{j}^{t}}\Big{|}<2^{t+1}\sigma,\quad
t=0,1,\cdots,2n-1,$ (3.12)
$\displaystyle\Big{|}{\sum_{j=1}^{n}\hat{a}_{j}\hat{g}_{j}^{t}-\sum_{j=1}^{n}a_{j}g_{j}^{t}}\Big{|}<2^{t+1}\sigma,\quad
t=0,1,\cdots,2n-1.$ (3.13)

Let

$\boldsymbol{\mathbf{d}}=\Big{(}\sum_{j=1}^{n}\hat{a}_{j}\hat{d}_{j}^{0}-\sum_{j=1}^{n}a_{j}d_{j}^{0},\quad\sum_{j=1}^{n}\hat{a}_{j}\hat{d}_{j}^{1}-\sum_{j=1}^{n}a_{j}d_{j}^{1},\quad\cdots,\quad\sum_{j=1}^{n}\hat{a}_{j}\hat{d}_{j}^{2n-1}-\sum_{j=1}^{n}a_{j}d_{j}^{2n-1}\Big{)}^{\top},$

and

$\boldsymbol{\mathbf{g}}=\Big{(}\sum_{j=1}^{n}\hat{a}_{j}\hat{g}_{j}^{0}-\sum_{j=1}^{n}a_{j}g_{j}^{0},\quad\sum_{j=1}^{n}\hat{a}_{j}\hat{g}_{j}^{1}-\sum_{j=1}^{n}a_{j}g_{j}^{1},\quad\cdots,\quad\sum_{j=1}^{n}\hat{a}_{j}\hat{g}_{j}^{2n-1}-\sum_{j=1}^{n}a_{j}g_{j}^{2n-1}\Big{)}^{\top}.$

Equations (3.12) and (3.13) imply respectively

$||\boldsymbol{\mathbf{d}}||_{2}<\frac{2^{2n+1}\sigma}{\sqrt{3}},\quad||\boldsymbol{\mathbf{g}}||_{2}<\frac{2^{2n+1}\sigma}{\sqrt{3}}.$

Note also that by (3.8), (3.9), (3.10), and (3.11), we get

$|\hat{d}_{j}|,|d_{j}|,|\hat{g}_{j}|,|g_{j}|\leq\sqrt{3},\quad j=1,\cdots,n.$

Define $d_{\min}:=\min_{p\neq q}\Big{|}{d_{p}-d_{q}}\Big{|}$ and
$g_{\min}:=\min_{p\neq q}\Big{|}{g_{p}-g_{q}}\Big{|}$. Applying Theorem 6.2,
we thus have that

$\Big{|}\Big{|}\eta_{n,n}(d_{1},\cdots,d_{n},\hat{d}_{1},\cdots,\hat{d}_{n})\Big{|}\Big{|}_{\infty}<\frac{(1+\sqrt{3})^{2n-1}}{d_{\min}^{n-1}}\frac{2^{2n+1}}{\sqrt{3}}\frac{\sigma}{m_{\min}},$
(3.14)

and

$\Big{|}\Big{|}\eta_{n,n}(g_{1},\cdots,g_{n},\hat{g}_{1},\cdots,\hat{g}_{n})\Big{|}\Big{|}_{\infty}<\frac{(1+\sqrt{3})^{2n-1}}{g_{\min}^{n-1}}\frac{2^{2n+1}}{\sqrt{3}}\frac{\sigma}{m_{\min}}.$
(3.15)

We now demonstrate that we can reorder $\hat{d}_{j},\hat{g}_{j}$ to have
$|\hat{d}_{j}-d_{j}|<\frac{d_{\min}}{2}$ and
$|\hat{g}_{j}-g_{j}|<\frac{g_{\min}}{2},j=1,\cdots,n$. First, since (3.8) and
(3.10) hold, by Lemma 3.2 we have

$\displaystyle d_{\min}\geq\frac{3}{2\pi}\min_{p\neq
q}\tau\Big{|}\Big{|}\boldsymbol{\mathbf{y}}_{p}-\boldsymbol{\mathbf{y}}_{q}\Big{|}\Big{|}_{1}\geq
11.475\Big{(}\frac{\sigma}{m_{\min}}\Big{)}^{\frac{1}{2n-1}}>2^{3/2}(1+\sqrt{3})\Big{(}\frac{2^{5/2}}{\sqrt{3}}\frac{\sigma}{m_{\min}}\Big{)}^{\frac{1}{2n-1}},$
(3.16)

and

$g_{\min}\geq\frac{3}{2\pi}\min_{p\neq
q}\tau\Big{|}\Big{|}\boldsymbol{\mathbf{y}}_{p}-\boldsymbol{\mathbf{y}}_{q}\Big{|}\Big{|}_{1}\geq
11.475\Big{(}\frac{\sigma}{m_{\min}}\Big{)}^{\frac{1}{2n-1}}>2^{3/2}(1+\sqrt{3})\Big{(}\frac{2^{5/2}}{\sqrt{3}}\frac{\sigma}{m_{\min}}\Big{)}^{\frac{1}{2n-1}},$

where we also use separation condition (2.5) in the above derivation. Let

$\epsilon_{d}=\frac{(1+\sqrt{3})^{2n-1}}{d_{\min}^{n-1}}\frac{2^{2n+1}}{\sqrt{3}}\frac{\sigma}{m_{\min}},\quad\epsilon_{g}=\frac{(1+\sqrt{3})^{2n-1}}{g_{\min}^{n-1}}\frac{2^{2n+1}}{\sqrt{3}}\frac{\sigma}{m_{\min}}.$

By (3.16), we have

$d_{\min}^{2n-1}\geq\frac{(1+\sqrt{3})^{2n-1}2^{3n+1}}{\sqrt{3}}\frac{\sigma}{m_{\min}},\
\text{or equivalently},d_{\min}^{n}\geq 2^{n}\epsilon_{d}.$

A similar result holds for $g_{\min}$ and $\epsilon_{g}$. Thus the conditions
of Lemma 6.8 are satisfied. By Lemma 6.8, we have that after reordering
$\hat{d}_{j},\hat{g}_{j}$,

$\Big{|}\hat{d}_{j}-d_{j}\Big{|}<\frac{d_{\min}}{2},\quad\Big{|}\hat{g}_{j}-g_{j}\Big{|}<\frac{g_{\min}}{2},$

and

$\displaystyle\Big{|}{\hat{d}_{j}-d_{j}}\Big{|}\leq\Big{(}\frac{2}{d_{\min}}\Big{)}^{n-1}\epsilon_{d}=\Big{(}\frac{1}{d_{\min}}\Big{)}^{2n-2}\frac{(1+\sqrt{3})^{2n-1}2^{3n}}{\sqrt{3}}\frac{\sigma}{m_{\min}},$
$\displaystyle\Big{|}{\hat{g}_{j}-g_{j}}\Big{|}\leq\Big{(}\frac{2}{g_{\min}}\Big{)}^{n-1}\epsilon_{g}=\Big{(}\frac{1}{g_{\min}}\Big{)}^{2n-2}\frac{(1+\sqrt{3})^{2n-1}2^{3n}}{\sqrt{3}}\frac{\sigma}{m_{\min}}.$

Observing

$e^{i\tau\boldsymbol{\mathbf{\hat{x}}}_{j,1}}-e^{i\tau\boldsymbol{\mathbf{x}}_{j,1}}=\frac{1}{2}\Big{(}\hat{d}_{j}-d_{j}+\hat{g}_{j}-g_{j}\Big{)},\qquad
e^{i\tau\boldsymbol{\mathbf{\hat{x}}}_{j,2}}-e^{i\tau\boldsymbol{\mathbf{x}}_{j,2}}=\frac{1}{2}\Big{(}\hat{d}_{j}-d_{j}-(\hat{g}_{j}-g_{j})\Big{)},$

we conclude that

$\Big{|}e^{i\tau\boldsymbol{\mathbf{\hat{x}}}_{j,1}}-e^{i\tau\boldsymbol{\mathbf{x}}_{j,1}}\Big{|}+\Big{|}e^{i\tau\boldsymbol{\mathbf{\hat{x}}}_{j,2}}-e^{i\tau\boldsymbol{\mathbf{x}}_{j,2}}\Big{|}\leq\Big{(}\Big{(}\frac{1}{d_{\min}}\Big{)}^{2n-2}+\Big{(}\frac{1}{g_{\min}}\Big{)}^{2n-2}\Big{)}\frac{(1+\sqrt{3})^{2n-1}2^{3n}}{\sqrt{3}}\frac{\sigma}{m_{\min}}.$

On the other hand, by (3.8) and (3.9),

$|\boldsymbol{\mathbf{\hat{x}}}_{j,1}-\boldsymbol{\mathbf{x}}_{j,1}|\leq\frac{\pi}{6}\quad\mbox{and}\quad|\boldsymbol{\mathbf{\hat{x}}}_{j,2}-\boldsymbol{\mathbf{x}}_{j,2}|\leq\frac{\pi}{6}.$

We further have

$\displaystyle\tau\Big{|}\boldsymbol{\mathbf{\hat{x}}}_{j,1}-\boldsymbol{\mathbf{x}}_{j,1}\Big{|}+\tau\Big{|}\boldsymbol{\mathbf{\hat{x}}}_{j,2}-\boldsymbol{\mathbf{x}}_{j,2}\Big{|}\leq\frac{\pi}{3}\Big{(}\Big{|}e^{i\boldsymbol{\mathbf{\hat{x}}}_{j,1}}-e^{i\boldsymbol{\mathbf{x}}_{j,1}}\Big{|}+\Big{|}e^{i\boldsymbol{\mathbf{\hat{x}}}_{j,2}}-e^{i\boldsymbol{\mathbf{x}}_{j,2}}\Big{|}\Big{)}$
$\displaystyle\leq$
$\displaystyle\Big{(}\Big{(}\frac{1}{d_{\min}}\Big{)}^{2n-2}+\Big{(}\frac{1}{g_{\min}}\Big{)}^{2n-2}\Big{)}\frac{(1+\sqrt{3})^{2n-1}2^{3n}\pi}{3\sqrt{3}}\frac{\sigma}{m_{\min}}.$

Recalling that $\tau=\frac{\Omega}{2n-1}$, we have

$\displaystyle\Big{|}\boldsymbol{\mathbf{\hat{x}}}_{j,1}-\boldsymbol{\mathbf{x}}_{j,1}\Big{|}+\Big{|}\boldsymbol{\mathbf{\hat{x}}}_{j,2}-\boldsymbol{\mathbf{x}}_{j,2}\Big{|}\leq\frac{2n-1}{\Omega}\Big{(}\Big{(}\frac{1}{d_{\min}}\Big{)}^{2n-2}+\Big{(}\frac{1}{g_{\min}}\Big{)}^{2n-2}\Big{)}\frac{(1+\sqrt{3})^{2n-1}2^{3n}\pi}{3\sqrt{3}}\frac{\sigma}{m_{\min}}.$

Note that by (3.16), we obtain that

$D_{\min}\leq\frac{2\pi(2n-1)}{3\Omega}d_{\min}\quad\mbox{and}\quad
D_{\min}\leq\frac{2\pi(2n-1)}{3\Omega}g_{\min}.$

Thus

$\displaystyle\Big{|}\Big{|}{\boldsymbol{\mathbf{\hat{x}}}_{j}-\boldsymbol{\mathbf{x}}_{j}}\Big{|}\Big{|}_{1}\leq$
$\displaystyle\frac{(1+\sqrt{3})^{2n-1}2^{3n+1}\pi(2n-1)}{3\sqrt{3}\Omega}\Big{(}\frac{2(2n-1)}{3}\Big{)}^{2n-2}\Big{(}\frac{\pi}{\Omega
D_{\min}}\Big{)}^{2n-2}\frac{\sigma}{m_{\min}}$ $\displaystyle=$
$\displaystyle\frac{(1+\sqrt{3})^{2n-1}2^{5n-1}(2n-1)^{2n-1}\pi}{3^{2n-0.5}\Omega}\Big{(}\frac{\pi}{\Omega
D_{\min}}\Big{)}^{2n-2}\frac{\sigma}{m_{\min}}.$

Since
$||\boldsymbol{\mathbf{\hat{y}}}_{j}-\boldsymbol{\mathbf{y}}_{j}||_{1}=||\boldsymbol{\mathbf{\hat{x}}}_{j}-\boldsymbol{\mathbf{x}}_{j}||_{1}$,
we further get

$\displaystyle\Big{|}\Big{|}{\boldsymbol{\mathbf{\hat{y}}}_{j}-\boldsymbol{\mathbf{y}}_{j}}\Big{|}\Big{|}_{1}\leq\frac{(1+\sqrt{3})^{2n-1}2^{5n-1}(2n-1)^{2n-1}\pi}{3^{2n-0.5}\Omega}\Big{(}\frac{\pi}{\Omega
D_{\min}}\Big{)}^{2n-2}\frac{\sigma}{m_{\min}}.$

Since
$D_{\min}\geq\frac{15.3\pi(n-0.5)}{\Omega}\Big{(}\frac{\sigma}{m_{\min}}\Big{)}^{\frac{1}{2n-1}}$,
together with the above estimate, we can also show that

$\Big{|}\Big{|}{\boldsymbol{\mathbf{\hat{y}}}_{j}-\boldsymbol{\mathbf{y}}_{j}}\Big{|}\Big{|}_{1}<\frac{D_{\min}}{2}.$

This completes the proof. ∎

##  4 An algorithm for the model order detection in two-dimensional DOA
estimation

In this section, based on the observations made in Section 3.1, we propose a
new algorithm, named coordinate-combination-based sweeping singular-value-
thresholding number detection algorithm, for the model order detection in two-
dimensional DOA estimations.

Figure 4.1: The geometry of a uniform rectangular array.

###  4.1 Problem formulation

The existing two-dimensional DOA algorithms usually try to estimate the
azimuth and elevation angles $(\theta_{j},\phi_{j})$’s that are shown in
Figure 4.1. More precisely, we consider $n$ narrowband signals/sources
impinging on an $(\Omega+1)\times(\Omega+1)$ uniform rectangular array (URA)
with $(\Omega+1)^{2}$ well calibrated and identically polarized antenna
elements. The signal received by these antenna elements in a single snapshot
can be expressed by

$\boldsymbol{\mathbf{Y}}(\boldsymbol{\mathbf{\omega}})=\sum_{j=1}^{n}s_{j}p_{j}e^{jkd_{x}\boldsymbol{\mathbf{\omega}}_{1}\boldsymbol{\mathbf{y}}_{j,1}}e^{jkd_{y}\boldsymbol{\mathbf{\omega}}_{2}\boldsymbol{\mathbf{y}}_{j,2}}+\boldsymbol{\mathbf{W}}(\boldsymbol{\mathbf{\omega}}),\quad\boldsymbol{\mathbf{\omega}}\in[0,1,\cdots,\Omega]^{2},$
(4.1)

where $s_{j}$ is the $j$-th incident signal, $p_{j}$ is a complex constant
denoting the signal/antenna polarization mismatch, $k$ represents the
wavenumber of the carrier frequency, and $d_{x}$ and $d_{y}$ denote the
distance between adjacent antenna element along the $x$-axis and $y$-axis,
respectively. $\boldsymbol{\mathbf{y}}_{j,1}=\sin\phi_{j}\cos\theta_{j}$ is
the direction component of signal $s_{j}$ propagating along the $x$-axis and
$\boldsymbol{\mathbf{y}}_{j,2}=\sin\phi_{j}\sin\theta_{j}$ is the one
propagating along the $y$-axis. The $\phi_{j}$ and $\theta_{j}$ denote
respectively the elevation and azimuth angles of $s_{j}$.
$\boldsymbol{\mathbf{W}}(\boldsymbol{\mathbf{\omega}})$ is the additive noise,
which is usually assumed to be white Gaussian noise.

For convenience, we consider the following simplified form of (4.1):

$\mathbf{Y}(\boldsymbol{\mathbf{\omega}})=\sum_{j=1}^{n}a_{j}e^{i\boldsymbol{\mathbf{y}}_{j}^{\top}\boldsymbol{\mathbf{\omega}}}+\mathbf{W}(\boldsymbol{\mathbf{\omega}}),\qquad\
\boldsymbol{\mathbf{\omega}}\in[0,1,\cdots,\Omega]^{2},$ (4.2)

where $\mathbf{W}$ is the noise with
$||\mathbf{W}(\boldsymbol{\mathbf{\omega}})||_{\infty}<\sigma$ and $\sigma$
being the noise level. We aim to recover stably the number of the signals and
the $\boldsymbol{\mathbf{y}}_{j}$’s, by which the elevation and azimuth angles
are stably resolved. For a better exposition, we still consider a discrete
measure $\mu=\sum_{j=1}^{n}a_{j}\delta_{\boldsymbol{\mathbf{y}}_{j}}$ and
denote the $a_{j}\delta_{\boldsymbol{\mathbf{y}}_{j}}$’s as sources. The
measurement (4.2) can be viewed as the noisy Fourier data of the measure $\mu$
at some discrete points.

In this section and the next one, we shall propose new algorithms for
detecting the model order and recovering the supports of $\mu$ from the
measurement (4.2). Our number detection method is based on thresholding on a
Hankel matrix assembled by data from modifications of (4.2) . The following
subsection shall introduce the details of the Hankel matrix formulation. We
refer to [2, 3, 74, 60, 56, 75, 36, 13, 28, 27, 43, 42] for other model
detecting algorithms.

###  4.2 Hankel matrix construction

The Hankel matrix is constructed by the following three steps.

Measurement modification by source translation
Due to the translation invariance, suppose the sources are supported in
$[-\lambda,\lambda]^{2}$, we consider them displacing with a vector
$\boldsymbol{\mathbf{v}}$ and get that
$\boldsymbol{\mathbf{x}}_{j}=\boldsymbol{\mathbf{y}}_{j}+\boldsymbol{\mathbf{v}}$.
Using a simple measurement modification technique, we obtain the measurement
for the new source
$\tilde{\mu}=\sum_{j=1}^{n}a_{j}\delta_{\boldsymbol{\mathbf{x}}_{j}}$.
Specifically, we consider

$\displaystyle\boldsymbol{\mathbf{X}}(\boldsymbol{\mathbf{\omega}})=$
$\displaystyle
e^{i\boldsymbol{\mathbf{v}}^{\top}\boldsymbol{\mathbf{\omega}}}\boldsymbol{\mathbf{Y}}(\boldsymbol{\mathbf{\omega}})=\sum_{j=1}^{n}a_{j}e^{i(\boldsymbol{\mathbf{y}}_{j}+\boldsymbol{\mathbf{v}})^{\top}\boldsymbol{\mathbf{\omega}}}+e^{i\boldsymbol{\mathbf{v}}^{\top}\boldsymbol{\mathbf{\omega}}}\mathbf{W}(\boldsymbol{\mathbf{\omega}})$
(4.3) $\displaystyle=$
$\displaystyle\sum_{j=1}^{n}a_{j}e^{i\boldsymbol{\mathbf{x}}_{j}^{\top}\boldsymbol{\mathbf{\omega}}}+\boldsymbol{\mathbf{\tilde{W}}}(\boldsymbol{\mathbf{\omega}}),\quad\boldsymbol{\mathbf{\omega}}\in[0,1,\cdots,\Omega]^{2},$

with $|\mathbf{\tilde{W}}(\boldsymbol{\mathbf{\omega}})|<\sigma$.

Measurement modification by coordinate-combination
The second procedure consists in modifying the measurement based on
coordinate-combination. For $s>0$, let $r=\frac{\Omega}{2s}$. From the
measurement $\boldsymbol{\mathbf{X}}$, we construct a list of new data given
by

$\boldsymbol{\mathbf{D}}(t)=\sum_{t_{1}+t_{2}=t,0\leq t_{1},t_{2}\leq
t}{t\choose
t_{1}}\boldsymbol{\mathbf{X}}(\boldsymbol{\mathbf{\omega}}_{rt_{1},rt_{2}}),\quad
t=0,\cdots,2s,$

where $\boldsymbol{\mathbf{\omega}}_{rt_{1},rt_{2}}=(rt_{1},rt_{2})^{\top}$.
Note that

$\displaystyle\boldsymbol{\mathbf{D}}(t)=$
$\displaystyle\sum_{j=1}^{n}a_{j}(e^{i\boldsymbol{\mathbf{x}}_{j,1}r}+e^{i\boldsymbol{\mathbf{x}}_{j,2}r})^{t}+\sum_{t_{1}+t_{2}=t,0\leq
t_{1},t_{2}\leq t}{t\choose
t_{1}}\boldsymbol{\mathbf{\tilde{W}}}(\boldsymbol{\mathbf{\omega}}_{rt_{1},rt_{2}})$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{n}a_{j}(e^{i\boldsymbol{\mathbf{x}}_{j,1}r}+e^{i\boldsymbol{\mathbf{x}}_{j,2}r})^{t}+\boldsymbol{\mathbf{\hat{W}}}(t),$

where $\boldsymbol{\mathbf{\hat{W}}}(t)=\sum_{t_{1}+t_{2}=t,0\leq
t_{1},t_{2}\leq t}{t\choose
t_{1}}\boldsymbol{\mathbf{\tilde{W}}}(\boldsymbol{\mathbf{\omega}}_{rt_{1},rt_{2}})$.

Hankel matrix construction and singular value decomposition
Finally, from these $\boldsymbol{\mathbf{D}}(t)$’s, we assemble the following
Hankel matrix

$\mathbf{H}(s)=\begin{pmatrix}\mathbf{D}(0)&\mathbf{D}(1)&\cdots&\mathbf{D}(s)\\\
\mathbf{D}(1)&\mathbf{D}(2)&\cdots&\mathbf{D}(s+1)\\\
\cdots&\cdots&\ddots&\cdots\\\
\mathbf{D}(s)&\mathbf{D}(s+1)&\cdots&\mathbf{D}(2s)\end{pmatrix}.$ (4.4)

We observe that $\mathbf{H}(s)$ has the decomposition

$\mathbf{H}(s)=BAB^{T}+\Delta,$ (4.5)

where $A=\text{diag}(a_{1},\cdots,a_{n})$ and
$B=\big{(}\phi_{s}(e^{i\boldsymbol{\mathbf{x}}_{j,1}r}+e^{i\boldsymbol{\mathbf{x}}_{j,2}r}),\cdots,\phi_{s}(e^{i\boldsymbol{\mathbf{x}}_{j,1}r}+e^{i\boldsymbol{\mathbf{x}}_{j,2}r})\big{)}$
with $\phi_{s}(\omega)$ being defined as

$\phi_{s}(\omega)=(1,\omega,\cdots,\omega^{s})^{\top},$ (4.6)

and

$\Delta=\begin{pmatrix}\mathbf{\hat{W}}(0)&\mathbf{\hat{W}}(1)&\cdots&\mathbf{\hat{W}}(s)\\\
\mathbf{\hat{W}}(1)&\mathbf{\hat{W}}(2)&\cdots&\mathbf{\hat{W}}(s+1)\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathbf{\hat{W}}(s)&\mathbf{\hat{W}}(s+1)&\cdots&\mathbf{\hat{W}}(2s)\end{pmatrix}.$
(4.7)

We denote the singular value decomposition of $\mathbf{H}(s)$ as

$\mathbf{H}(s)=\hat{U}\hat{\Sigma}\hat{U}^{*},$

where
$\hat{\Sigma}=\text{diag}(\hat{\sigma}_{1},\cdots,\hat{\sigma}_{n},\hat{\sigma}_{n+1},\cdots,\hat{\sigma}_{s+1})$
with the singular values $\hat{\sigma}_{j}$, $1\leq j\leq s+1$, ordered in a
decreasing manner. The source number $n$ is then detected by a thresholding on
these singular values. In the next subsection we will provide the theoretical
guarantee of the threshold.

###  4.3 Theoretical guarantee

Note that when there is no noise, $\mathbf{H}(s)=BAB^{\top}$. We have the
following estimate for the singular values of $BAB^{\top}$.

######  Lemma 4.1.

Let $n\geq 2$, $s\geq n$,
$\boldsymbol{\mathbf{y}}_{j}\in[-\frac{s\pi}{6\Omega},\frac{s\pi}{6\Omega}]^{2},1\leq
j\leq n,$ and $\boldsymbol{\mathbf{v}}$ in (4.3) be
$(0,\frac{s\pi}{\Omega})^{\top}$. Let

$\sigma_{1},\cdots,\sigma_{n},0,\cdots,0$

be the singular values of $BAB^{T}$ in (4.5) ordered in a decreasing manner.
Then the following estimate holds

$\sigma_{n}\geq\frac{m_{\min}\big{(}3\theta_{\min}(\Omega,s)\big{)}^{2n-2}}{n(2(1+\sqrt{3})\pi)^{2n-2}},$
(4.8)

where $\theta_{\min}(\Omega,s)=\min_{p\neq
j}\Big{|}\Big{|}{\boldsymbol{\mathbf{y}}_{p}\frac{\Omega}{2s}-\boldsymbol{\mathbf{y}}_{j}\frac{\Omega}{2s}}\Big{|}\Big{|}_{1}$.

###### Proof.

Recall that $\sigma_{n}$ is the minimum nonzero singular value of
$BAB^{\top}$. Let $\ker(B^{\top})$ be the kernel space of $B^{\top}$ and
$\ker^{\perp}(B^{\top})$ be its orthogonal complement. Then we have

$\displaystyle\sigma_{n}=\min_{||x||_{2}=1,x\in\ker^{\perp}(B^{\top})}||BAB^{\top}x||_{2}\geq\sigma_{\min}(BA)\sigma_{n}(B^{\top})$
$\displaystyle\geq$
$\displaystyle\sigma_{\min}(B)\sigma_{\min}(A)\sigma_{\min}(B).$

On the other hand, since by the condition of the lemma
$\boldsymbol{\mathbf{x}}_{j}=\boldsymbol{\mathbf{y}}_{j}+\boldsymbol{\mathbf{v}}\in[-\frac{s\pi}{6\Omega},\frac{s\pi}{6\Omega}]\times[\frac{5s\pi}{6\Omega},\frac{7s\pi}{6\Omega}]$,
we have
$\frac{\Omega\boldsymbol{\mathbf{x}}_{j}}{2s}\in[-\frac{\pi}{12},\frac{\pi}{12}]\times[\frac{5\pi}{12},\frac{7\pi}{12}]$.
Thus, by Lemma 3.2, for $r=\frac{\Omega}{2s}$,

$\min_{p\neq
q}\Big{|}{e^{i\boldsymbol{\mathbf{x}}_{p,1}r}+e^{i\boldsymbol{\mathbf{x}}_{p,2}r}-(e^{i\boldsymbol{\mathbf{x}}_{q,1}r}+e^{i\boldsymbol{\mathbf{x}}_{q,2}r})}\Big{|}\geq\frac{3}{2\pi}\theta_{\min}(\Omega,s).$

Note also that
$|e^{i\boldsymbol{\mathbf{x}}_{p,1}r}+e^{i\boldsymbol{\mathbf{x}}_{p,2}r}|\leq\sqrt{3}$.
Thus applying Lemma 6.3 and Corollary 6.2, we have

$\displaystyle\sigma_{\min}(B)\geq\frac{1}{\sqrt{n}}\frac{\big{(}\frac{3}{2\pi}\theta_{\min}(\Omega,s)\big{)}^{n-1}}{(1+\sqrt{3})^{n-1}}.$

Then, it follows that

$\displaystyle\sigma_{n}\geq\sigma_{\min}(A)\Big{(}\frac{\big{(}\frac{3}{2\pi}\theta_{\min}(\Omega,s)\big{)}^{n-1}}{(1+\sqrt{3})^{n-1}}\Big{)}^{2}\geq\frac{m_{\min}\big{(}3\theta_{\min}(\Omega,s)\big{)}^{2n-2}}{n(2(1+\sqrt{3})\pi)^{2n-2}}.$

∎

We now present the main result on the threshold for the singular values of the
matrix $\mathbf{H}(s)$.

######  Theorem 4.1.

Let $n\geq 2,s\geq n$ and
$\mu=\sum_{j=1}^{n}a_{j}\delta_{\boldsymbol{\mathbf{y}}_{j}}$ with
$\boldsymbol{\mathbf{y}}_{j}\in[-\frac{s\pi}{6\Omega},\frac{s\pi}{6\Omega}]^{2},1\leq
j\leq n$. Let $\boldsymbol{\mathbf{v}}$ in (4.3) be equal to
$(0,\frac{s\pi}{\Omega})^{\top}$. Then for the singular values of
$\boldsymbol{\mathbf{H}}(s)$ in (4.4), We have

$\hat{\sigma}_{j}<\frac{4^{s+1}\sigma}{3},\quad j=n+1,\cdots,s+1.$ (4.9)

Moreover, if the following separation condition is satisfied

$\min_{p\neq
j}\Big{|}\Big{|}{\boldsymbol{\mathbf{y}}_{p}-\boldsymbol{\mathbf{y}}_{j}}\Big{|}\Big{|}_{1}\geq\frac{4(1+\sqrt{3})\pi
s}{3\Omega}\Big{(}\frac{2n4^{s+1}}{3}\frac{\sigma}{m_{\min}}\Big{)}^{\frac{1}{2n-2}},$
(4.10)

then

$\hat{\sigma}_{n}\geq\frac{4^{s+1}\sigma}{3}.$ (4.11)

###### Proof.

We first estimate $||\Delta||_{2}$ for $\Delta$ in (4.7). By the definition of
$\boldsymbol{\mathbf{\hat{W}}}(t)$, we have
$|\boldsymbol{\mathbf{\hat{W}}}(t)<2^{t}\sigma$. Thus
$||\Delta||_{2}\leq||\Delta||_{F}<\frac{4^{s+1}\sigma}{3}$. By Weyl’s theorem,
we have $|\hat{\sigma}_{j}-\sigma_{j}|\leq||\Delta||_{2},j=1,\cdots,n$.
Together with $\sigma_{j}=0,n+1\leq j\leq s+1$, we get
$|\hat{\sigma}_{j}|\leq||\Delta||_{2}<\frac{4^{s+1}\sigma}{3},n+1\leq j\leq
s+1$. This proves (4.9).

Let $\theta_{\min}(\Omega,s)=\frac{\Omega}{2s}\min_{p\neq
q}\Big{|}\Big{|}{\boldsymbol{\mathbf{y}}_{p}-\boldsymbol{\mathbf{y}}_{q}}\Big{|}\Big{|}_{1}$.
The separation condition (4.10) implies that

$\theta_{\min}(\Omega,s)\geq\frac{2(1+\sqrt{3})\pi}{3}\Big{(}\frac{2n4^{s+1}}{3}\frac{\sigma}{m_{\min}}\Big{)}^{\frac{1}{2n-2}}.$

By Lemma 4.1, we have

$\displaystyle\sigma_{n}\geq\frac{m_{\min}\big{(}3\theta_{\min}(\Omega,s)\big{)}^{2n-2}}{n(2(1+\sqrt{3})\pi)^{2n-2}}>2\frac{4^{s+1}\sigma}{3}.$
(4.12)

Similarly, by Weyl’s theorem,
$|\hat{\sigma}_{n}-\sigma_{n}|\leq||\Delta||_{2}$. Thus, $\hat{\sigma}_{n}\geq
2(s+1)\sigma-||\Delta||_{2}\geq\frac{4^{s+1}\sigma}{3}$. The conclusion (4.11)
then follows. ∎

###  4.4 Coordinate-combination-based sweeping singular-value-thresholding
number detection algorithm

Based on Theorem 4.1, we can propose a simple thresholding algorithm,
Algorithm 1, for the number detection.

Input: Number $s$; Noise level $\sigma$;

Input: Measurement:
$\mathbf{Y}(\boldsymbol{\mathbf{\omega}}),\boldsymbol{\mathbf{\omega}}\in[0,1,\cdots,\Omega]^{2}$;

Input: Translation vector $\boldsymbol{\mathbf{v}}$ in $\mathbb{R}^{2}$;

1: Modify the measurement and get
$\boldsymbol{\mathbf{X}}(\boldsymbol{\mathbf{\omega}})=e^{i\boldsymbol{\mathbf{v}}^{\top}\boldsymbol{\mathbf{\omega}}}\mathbf{Y}(\boldsymbol{\mathbf{\omega}})$;

2: Let $r=\Omega\mod 2s$, formulate
$\boldsymbol{\mathbf{D}}(t)=\sum_{t_{1}+t_{2}=t,0\leq t_{1},t_{2}\leq
t}{t\choose
t_{1}}\boldsymbol{\mathbf{X}}(\boldsymbol{\mathbf{\omega}}_{rt_{1},rt_{2}}),\quad
t=0,\cdots,2s$;

3: Assemble the $(s+1)\times(s+1)$ Hankel matrix $\mathbf{H}(s)$ like (4.4)
from $\boldsymbol{\mathbf{D}}(t)$’s, and compute the singular value of
$\mathbf{H}(s)$ as $\hat{\sigma}_{1},\cdots,\hat{\sigma}_{s+1}$ distributed in
a decreasing manner;

4: Determine $n$ by $\hat{\sigma}_{n}\geq\frac{4^{s+1}\sigma}{3}$ and
$\hat{\sigma}_{j}<\frac{4^{s+1}\sigma}{3},j=n+1,\cdots,s+1$;

Return: $n$

Algorithm 1 Coordinate-combination-based singular-value-thresholding number
detection algorithm

Note that for Algorithm 1 to work, in addition to the smallness of the noise
level $\sigma$, we also need the integer $s$ to be larger than the source
number. However, a suitable $s$ is not easy to estimate and large $s$ may
incur a deterioration of the resolution as indicated by (4.10). To remedy this
issue, we propose a sweeping singular-value-thresholding number detection
algorithm (Algorithm 2) below. In short, we detect the number $n_{recover}$ by
Algorithm 1 for all $s$ from $2$ to $\lfloor\frac{\Omega-1}{2}\rfloor$, and
choose the greatest one $n_{\max}$ as the number of point sources. When the
detected $n_{recover}$ becomes smaller than $n_{\max}$ for a large number of
iterations, we will stop the loop. The details are summarized in Algorithm2
below.

We remark that when $s=n$ and the point sources satisfy

$\min_{p\neq
q}\Big{|}\Big{|}{\boldsymbol{\mathbf{y}}_{p}-\boldsymbol{\mathbf{y}}_{q}}\Big{|}\Big{|}_{1}\geq\frac{Cn}{\Omega}\Big{(}\frac{\sigma}{m_{\min}}\Big{)}^{\frac{1}{2n-2}},$
(4.13)

for some constant $C$, then (4.10) is satisfied. Thus by Theorem 4.1, for a
suitable choice of $\boldsymbol{\mathbf{v}}$, Algorithm 1 can exactly detect
the number $n$ when $s=n$. As $s$ increases to values greater than $n$, (4.9)
implies that the number detected by Algorithm 1 will not exceed $n$.
Therefore, the sweeping singular-value-thresholding algorithm (Algorithm 2)
can detect the exact number $n$ when $\Omega$ is greater than $2n+1$ and the
point sources are separated by the minimal separation distance we derived in
Theorem 2.1. This demonstrates the optimal performance of Algorithm 2. We also
remark that the theoretical threshold derived in Theorem 4.1 seems to be
larger than the one that is needed. One can improve the algorithm by choosing
smaller threshold. Deriving new estimates for the thresholds in different
cases is another interesting problem.

Input: Noise level $\sigma$; Measurement: $\mathbf{Y}$; Translation vector
$\boldsymbol{\mathbf{v}}$;

Input: $n_{max}=0$, $smax_{index}=2$

for  _$s=2:\lfloor\frac{\Omega-1}{2}\rfloor$_ do

Input $s,\sigma,\mathbf{Y},\boldsymbol{\mathbf{v}}$ to Algorithm 1, save the
output of Algorithm 1 as $n_{recover}$;

if  _$n_{recover} >n_{max}$_ then

$n_{max}=n_{recover}$;

$smax_{index}=s$;

if  _$s\geq$ $smax_{index}+2$_ then

break;

Return $n_{max}$.

Algorithm 2 Coordinate-combination-based sweeping singular-value-thresholding
number detection algorithm

###  4.5 Phase transition and performance of Algorithm 2

In this subsection, we conduct numerical experiments to demonstrate the phase
transition phenomenon regarding the super-resolution factor (SRF) and the SNR
using Algorithm 2. We consider recovering the number of three and four
sources. We fix $\Omega=10$ and detect the source number from their noisy
Fourier data at $[0,1,\cdots,\Omega]^{2}$. We consider sources in
$[0,\frac{\pi}{2}]^{2}$ and the translation vector in Algorithm 1 is
$\boldsymbol{\mathbf{v}}=(0,\frac{\pi}{2})^{\top}$. The noise level is
$\sigma$ and the minimum separation distance between sources is $D_{\min}$. We
perform $10000$ random experiments (the randomness is in the choice of
($d_{\min}$,$\sigma$, $\boldsymbol{\mathbf{y}}_{j}$, $a_{j}$)) and detect the
source number by Algorithm 2. We record the number of each successful
detection (source number is detected exactly) and failed detection. Figures
4.2 shows the result for the successful and unsuccessfully recovery in the
parameter space $\log(SNR)$ versus $\log(SRF)$ . It is observed that there is
a line with slope ($2n-2$) in the parameter space of $\log(SRF)$-$\log(SNR)$
above which the number detection is always successful. This phase transition
phenomenon is exactly the one predicted by our theoretical results in Theorems
2.1 and 4.1. It also illustrates the efficiency of Algorithm 2 as it can
resolve the source number correctly in the regime where the source separation
distance is of the order of the computational resolution limit.

(a) detection success

(b) detection success

Figure 4.2: Plots of the successful and the unsuccessful number detection by
Algorithm 2 depending on the relation between $\log(SRF)$ and
$\log(\frac{1}{\sigma})$. (a) illustrates that three sources can be exactly
detected if $\log(\frac{1}{\sigma})$ is above a line of slope $4$ in the
parameter space. (b) illustrates that four sources can be exactly detected if
$\log(\frac{1}{\sigma})$ is above a line of slope $6$ in the parameter space.

##  5 An algorithm for the source reconstruction in two-dimensional DOA
problems

In this section, based on the idea of coordinate-combination, we propose a new
MUSIC algorithm for resolving the sources in the two-dimensional DOA
estimation. Our algorithm is named as coordinate-combination-based MUSIC
algorithm; see Algorithm 4.

###  5.1 Hankel matrix construction

Similarly to the number detection algorithm in the above section, the MUSIC
algorithm also relies on a singular value decomposition of certain Hankel
matrix. Compared to conventional MUSIC-based DOA algorithms, the main novelty
of our algorithm lies in a different way of assembling Hankel matrices.
Similarly to Section 4.2, the Hankel matrix construction here is also based on
observations in Section 3.1 and the details are presented below.

Measurement modification by source translation
We consider the same model setting as (4.2) for the available measurement. We
also perform the source translation and modify the measurement to get

$\displaystyle\boldsymbol{\mathbf{X}}(\boldsymbol{\mathbf{\omega}})=$
$\displaystyle
e^{i\boldsymbol{\mathbf{v}}^{\top}\boldsymbol{\mathbf{\omega}}}\boldsymbol{\mathbf{Y}}(\boldsymbol{\mathbf{\omega}})=\sum_{j=1}^{n}a_{j}e^{i(\boldsymbol{\mathbf{y}}_{j}+\boldsymbol{\mathbf{v}})^{\top}\boldsymbol{\mathbf{\omega}}}+e^{i\boldsymbol{\mathbf{v}}^{\top}\boldsymbol{\mathbf{\omega}}}\mathbf{W}(\boldsymbol{\mathbf{\omega}})$
(5.1) $\displaystyle=$
$\displaystyle\sum_{j=1}^{n}a_{j}e^{i\boldsymbol{\mathbf{x}}_{j}^{\top}\boldsymbol{\mathbf{\omega}}}+\boldsymbol{\mathbf{\tilde{W}}}(\boldsymbol{\mathbf{\omega}}),\quad\boldsymbol{\mathbf{\omega}}\in[0,1,\cdots,\Omega]^{2},$

where
$\boldsymbol{\mathbf{x}}_{j}=\boldsymbol{\mathbf{y}}_{j}+\boldsymbol{\mathbf{v}}$
for a suitable $\boldsymbol{\mathbf{v}}\in\mathbb{R}^{2}$ and
$|\mathbf{\tilde{W}}(\boldsymbol{\mathbf{\omega}})|<\sigma$.

Measurement modification by the coordinate-combination technique
Let $s=\lfloor\frac{\Omega}{2}\rfloor$. From the modified measurement
$\boldsymbol{\mathbf{X}}(\boldsymbol{\mathbf{\omega}})$, we construct the
following two lists of data:

$\displaystyle\boldsymbol{\mathbf{D}}(t)=\sum_{t_{1}+t_{2}=t,0\leq
t_{1},t_{2}\leq t}{t\choose
t_{1}}\boldsymbol{\mathbf{X}}(\boldsymbol{\mathbf{\omega}}_{t_{1},t_{2}}),\quad
t=0,\cdots,2s,$
$\displaystyle\boldsymbol{\mathbf{G}}(t)=\sum_{t_{1}+t_{2}=t,0\leq
t_{1},t_{2}\leq t}(-1)^{t_{2}}{t\choose
t_{1}}\boldsymbol{\mathbf{X}}(\boldsymbol{\mathbf{\omega}}_{t_{1},t_{2}}),\quad
t=0,\cdots,2s,$

where $\boldsymbol{\mathbf{\omega}}_{t_{1},t_{2}}=(t_{1},t_{2})^{\top}$. Note
that

$\displaystyle\boldsymbol{\mathbf{D}}(t)=$
$\displaystyle\sum_{j=1}^{n}a_{j}(e^{i\boldsymbol{\mathbf{x}}_{j,1}}+e^{i\boldsymbol{\mathbf{x}}_{j,2}})^{t}+\sum_{t_{1}+t_{2}=t,0\leq
t_{1},t_{2}\leq t}{t\choose
t_{1}}\boldsymbol{\mathbf{\tilde{W}}}(\boldsymbol{\mathbf{\omega}}_{t_{1},t_{2}})$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{n}a_{j}(e^{i\boldsymbol{\mathbf{x}}_{j,1}}+e^{i\boldsymbol{\mathbf{x}}_{j,2}})^{t}+\boldsymbol{\mathbf{\hat{W}}}_{d}(t),$
$\displaystyle\boldsymbol{\mathbf{G}}(t)=$
$\displaystyle\sum_{j=1}^{n}a_{j}(e^{i\boldsymbol{\mathbf{x}}_{j,1}}-e^{i\boldsymbol{\mathbf{x}}_{j,2}})^{t}+\sum_{t_{1}+t_{2}=t,0\leq
t_{1},t_{2}\leq t}(-1)^{t_{2}}{t\choose
t_{1}}\boldsymbol{\mathbf{\tilde{W}}}(\boldsymbol{\mathbf{\omega}}_{t_{1},t_{2}})$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{n}a_{j}(e^{i\boldsymbol{\mathbf{x}}_{j,1}}-e^{i\boldsymbol{\mathbf{x}}_{j,2}})^{t}+\boldsymbol{\mathbf{\hat{W}}}_{g}(t),$

where $\boldsymbol{\mathbf{\hat{W}}}_{d}(t)=\sum_{t_{1}+t_{2}=t,0\leq
t_{1},t_{2}\leq t}{t\choose
t_{1}}\boldsymbol{\mathbf{\tilde{W}}}(\boldsymbol{\mathbf{\omega}}_{t_{1},t_{2}})$
and $\boldsymbol{\mathbf{\hat{W}}}_{g}(t)=\sum_{t_{1}+t_{2}=t,0\leq
t_{1},t_{2}\leq t}(-1)^{t_{2}}{t\choose
t_{1}}\boldsymbol{\mathbf{\tilde{W}}}(\boldsymbol{\mathbf{\omega}}_{t_{1},t_{2}})$.

Hankel matrix construction
Finally, from these $\boldsymbol{\mathbf{D}}(t),\boldsymbol{\mathbf{G}}(t)$’s,
we assemble the following Hankel matrices:

$\mathbf{H}_{d}(s)=\begin{pmatrix}\mathbf{D}(0)&\mathbf{D}(1)&\cdots&\mathbf{D}(s)\\\
\mathbf{D}(1)&\mathbf{D}(2)&\cdots&\mathbf{D}(s+1)\\\
\cdots&\cdots&\ddots&\cdots\\\
\mathbf{D}(s)&\mathbf{D}(s+1)&\cdots&\mathbf{D}(2s)\end{pmatrix},\quad\mathbf{H}_{g}(s)=\begin{pmatrix}\mathbf{G}(0)&\mathbf{G}(1)&\cdots&\mathbf{G}(s)\\\
\mathbf{G}(1)&\mathbf{G}(2)&\cdots&\mathbf{G}(s+1)\\\
\cdots&\cdots&\ddots&\cdots\\\
\mathbf{G}(s)&\mathbf{G}(s+1)&\cdots&\mathbf{G}(2s)\end{pmatrix}.$ (5.2)

###  5.2 Standard MUSIC algorithm

In this subsection, we perform the standard MUSIC algorithm [58, 65, 40, 45]
for the Hankel matrix
$\boldsymbol{\mathbf{H}}_{d}(s),\boldsymbol{\mathbf{H}}_{g}(s)$ in (5.2). For
ease of presentation, we only introduce the MUSIC algorithm for
$\boldsymbol{\mathbf{H}}_{d}(s)$. The one for $\boldsymbol{\mathbf{H}}_{g}(s)$
can be developed in the same manner. Our algorithm first performs the singular
value decomposition of $\boldsymbol{\mathbf{H}}_{d}(s)$,

$\boldsymbol{\mathbf{H}}_{d}(s)=\hat{U}\hat{\Sigma}\hat{U}^{*}=[\hat{U}_{1}\quad\hat{U}_{2}]\text{diag}(\hat{\sigma}_{1},\hat{\sigma}_{2},\cdots,\hat{\sigma}_{n},\hat{\sigma}_{n+1},\cdots,\hat{\sigma}_{s+1})[\hat{U}_{1}\quad\hat{U}_{2}]^{*},$

where
$\hat{U}_{1}=(\hat{U}(1),\cdots,\hat{U}(n)),\hat{U}_{2}=(\hat{U}(n+1),\cdots,\hat{U}(s+1))$
with $n$ being the estimated source number (model order). The source number
$n$ can be detected by Algorithm 2 and many other algorithms such as those in
[2, 74, 60, 75, 13, 28, 27, 43, 42]. Denote the orthogonal projection onto the
space $\hat{U}_{2}$ by $\hat{P}_{2}x=\hat{U}_{2}(\hat{U}_{2}^{*}x)$. For a
test vector $\Phi(d)=(1,d,\cdots,d^{s})^{\top}$, one defines the MUSIC imaging
functional

$\displaystyle\hat{J}(d)=\frac{||\Phi(d)||_{2}}{||\hat{P}_{2}\Phi(d)||_{2}}=\frac{||\Phi(d)||_{2}}{||\hat{U}_{2}^{*}\Phi(d)||_{2}}.$

The local maximizer of $\hat{J}(d)$ indicates the supports of the sources. In
practice, one can test evenly spaced points in a specified region and plot the
discrete imaging functional and then determine the sources by detecting the
peaks. In our case, we only need to test some discrete points $d\in\mathbb{C}$
with $|d|\leq 2$ and select the peak by certain algorithms (such as the one in
[45] or its two-dimensional analog). Finally, we summarize the standard MUSIC
algorithm in Algorithm 3 below.

Input: Source number $n$;

Input: Modified measurements: $\boldsymbol{\mathbf{D}}(t)$ (or
$\boldsymbol{\mathbf{G}}(t)$), $t=0,\cdots,s$ with $s\geq n$;

Input: Test points $d$’s;

1: Formulate the $(s+1)\times(s+1)$ Hankel matrix
$\boldsymbol{\mathbf{H}}_{d}(s)$ from $\boldsymbol{\mathbf{D}}(t)$’s as (5.2);

2: Compute the singular vectors of $\boldsymbol{\mathbf{H}}_{d}(s)$ as
$\hat{U}(1),\hat{U}(2),\cdots,\hat{U}(s+1)$ and form the noise space
$\hat{U}_{2}=(\hat{U}(n+1),\cdots,\hat{U}(s+1))$;

3: For test points $d$’s, construct the test vector
$\Phi(d)=(1,d,\cdots,d^{s})^{\top}$;

4: Plot the MUSIC imaging functional
$\hat{J}(d)=\frac{||\Phi(d)||_{2}}{||\hat{U}_{2}^{*}\Phi(d)||_{2}}$;

5: Select the peak locations $\hat{d}_{j}$’s in the plot of $\hat{J}(d)$.

Algorithm 3 Standard MUSIC algorithm

###  5.3 Coordinate-combination-based MUSIC algorithm

After applying the MUSIC algorithm to both
$\boldsymbol{\mathbf{H}}_{d}(s),\boldsymbol{\mathbf{H}}_{g}(s)$, we expect to
reconstruct $n$ $\hat{d}_{j}$’s which is close to
$d_{j}=e^{i\boldsymbol{\mathbf{x}}_{j,1}}+e^{i\boldsymbol{\mathbf{x}}_{j,2}}$,
and $n$ $\hat{g}_{j}$’s which is close to
$g_{j}=e^{i\boldsymbol{\mathbf{x}}_{j,1}}-e^{i\boldsymbol{\mathbf{x}}_{j,2}}$.
The next question is how to link the pair $\hat{d}_{j},\hat{g}_{j}$ that
correspond to the same source. This is an inevitable pair matching issues in
most of the two-dimensional DOA algorithms [42], where ad hoc schemes [79, 34,
14, 77] were derived to associate the estimated azimuth and elevation angles.
Here, in contrast with conventional DOA algorithms, we do not need to link the
azimuth and elevation angles but to link $\hat{d}_{j}$ and $\hat{g}_{j}$.

Observe that $|d_{j}+g_{j}|=|2e^{i\boldsymbol{\mathbf{x}}_{j,1}}|=2$ and
$|d_{j}-g_{j}|=|2e^{i\boldsymbol{\mathbf{x}}_{j,2}}|=2$. We can use this
criterion to match the pair $\hat{d}_{j},\hat{g}_{j}$ that they should satisfy

$|\hat{d}_{j}+\hat{g}_{j}|\approx 2,\quad|\hat{d}_{j}-\hat{g}_{j}|\approx 2.$
(5.3)

For example, we could consider the following minimization problem:

$\min_{\pi\in\zeta(n)}\sum_{j=1}^{n}\Big{|}{|\hat{d}_{j}+\hat{g}_{\pi_{j}}|-2}\Big{|}+\Big{|}{|\hat{d}_{j}-\hat{g}_{\pi_{j}}|-2}\Big{|},$
(5.4)

where $\zeta(n)$ is the set of all permutations of $\\{1,\cdots,n\\}$. This
can be viewed as a balanced assignment problem [52], which can be solved
efficiently by many algorithms such as the Hungarian algorithm.

We remark that our pair matching algorithm is not the one usually required in
other one-dimensional based DOA algorithms. Unlike our case, the other pair
matching problem is not an assignment problem, wherefore the pair matching is
usually time consuming or complex processing is conducted to reduce the
computational cost.

Input: Source number $n$; noise level $\sigma$;

Input: Measurement:
$\mathbf{Y}(\boldsymbol{\mathbf{\omega}}),\boldsymbol{\mathbf{\omega}}\in[0,1,\cdots,\Omega]^{2}$;

Input: Translation vector $\boldsymbol{\mathbf{v}}$ in $\mathbb{R}^{2}$;

Input: Evenly spaced test points $d\in\mathbb{C}$ with $|d|\leq 2$;

1: Modify the measurement and get
$\boldsymbol{\mathbf{X}}(\boldsymbol{\mathbf{\omega}})=e^{i\boldsymbol{\mathbf{v}}^{\top}\boldsymbol{\mathbf{\omega}}}\mathbf{Y}(\boldsymbol{\mathbf{\omega}})$;

2: Let $s=\lfloor\frac{\Omega}{2}\rfloor$, formulate
$\boldsymbol{\mathbf{D}}(t)=\sum_{t_{1}+t_{2}=t,0\leq t_{1},t_{2}\leq
t}{t\choose
t_{1}}\boldsymbol{\mathbf{X}}(\boldsymbol{\mathbf{\omega}}_{t_{1},t_{2}}),\quad\boldsymbol{\mathbf{G}}(t)=\sum_{t_{1}+t_{2}=t,0\leq
t_{1},t_{2}\leq t}(-1)^{t_{2}}{t\choose
t_{1}}\boldsymbol{\mathbf{X}}(\boldsymbol{\mathbf{\omega}}_{t_{1},t_{2}}),\quad
t=0,\cdots,2s$;

3: Input $\boldsymbol{\mathbf{D}},n$ and test points $d$’s into Algorithm 3
and get the output $\hat{d}_{1},\cdots,\hat{d}_{n}$;

4: Input $\boldsymbol{\mathbf{G}},n$ and test points $d$’s into Algorithm 3
and get the output $\hat{g}_{1},\cdots,\hat{g}_{n}$;

5: Matching the $\hat{d}_{j},\hat{g}_{j}$’s by applying an assignment
algorithm (match pairs in matlab) to solve (5.4) and get the pair list
$\\{(\hat{d}_{j},\hat{g}_{j})\\}^{j=1^{n}}$;

6: Get $\frac{\hat{d}_{j}+\hat{g}_{j}}{2}$ and
$\frac{\hat{d}_{j}-\hat{g}_{j}}{2},j=1,\cdots,n$. Get
$e^{i\boldsymbol{\mathbf{\hat{x}}}_{j,1}r}$ by considering the closest point
to $\frac{\hat{d}_{j}+\hat{g}_{j}}{2}$ on the unit circle. Get
$e^{i\boldsymbol{\mathbf{\hat{x}}}_{j,2}}$ by considering the closest point to
$\frac{\hat{d}_{j}-\hat{g}_{j}}{2}$ on the unit circle;

7: The recovered
$\boldsymbol{\mathbf{\hat{x}}}_{j}=(\boldsymbol{\mathbf{\hat{x}}}_{j,1},\boldsymbol{\mathbf{\hat{x}}}_{j,2})^{\top}$.
Reconstruct
$\boldsymbol{\mathbf{\hat{y}}}_{j}=\boldsymbol{\mathbf{\hat{x}}}_{j}-\boldsymbol{\mathbf{v}},j=1,\cdots,n$;

Return:
$\boldsymbol{\mathbf{\hat{y}}}_{1},\cdots,\boldsymbol{\mathbf{\hat{y}}}_{n}$.

Algorithm 4 Coordinate-combination-based MUSIC algorithm for two-dimensional
DOA

###  5.4 Superiority of the algorithm

####  5.4.1 Overcome the issue of separation distance loss in conventional
two-dimensional DOA algorithms

Despite the fact that different recovering methods are proposed for DOA
estimation in two dimensions, the conventional way for tackling the problem
has hardly exceeded the scope of recovering the two direction (x- and
y-direction) components of sources individually. Thus, as illustrated in
Figure 5.1, severe loss of the source separation distance in one dimension is
always an inevitable issue that causes unstable recovery of the direction
components. Most of the researches ignored this issue and some papers [73, 72]
proposed ad hoc schemes to enhance the reconstruction but in a complex manner.

Our method is a new one-dimensional-based algorithm where the issue of severe
source separation distance loss is avoided in a simple way. In our algorithm,
the separation distance between direction components of sources are still
preserved. This has been demonstrated by Lemma 3.2 for
$\boldsymbol{\mathbf{\theta}}_{j}\in\mathbb{[}0,\pi]^{2},j=1,2$ with
$\frac{\pi}{3}\leq\boldsymbol{\mathbf{\theta}}_{j,2}-\boldsymbol{\mathbf{\theta}}_{j,1}\leq\frac{2}{3}\pi,j=1,2$.
Furthermore, Theorem 5.1 shows that, for
$\boldsymbol{\mathbf{y}}_{j}\in[0,\frac{\pi}{2}]^{2}$ and
$\boldsymbol{\mathbf{v}}=(0,\frac{\pi}{2})^{\top}$, the separation distance
between
$\boldsymbol{\mathbf{x}}_{j}=\boldsymbol{\mathbf{y}}_{j}+\boldsymbol{\mathbf{v}}$’s
can be preserved after the coordinate-combination. By Theorem 5.1, if the
distance between the $\boldsymbol{\mathbf{x}}_{j}$’s is a certain constant
$C$, then the distance between
$e^{i\boldsymbol{\mathbf{x}}_{j,1}}+e^{i\boldsymbol{\mathbf{x}}_{j,2}}$ is
larger than $\frac{2C}{\pi^{2}}$ times the original distance. For better
results of preservation of the distance, as indicated by Theorems 2.1 and 2.2,
we could consider sources in a smaller region with a specified translation. In
the numerical experiments presented in this paper, for ease of discussion and
presentation, we will consider sources in $[0,\frac{\pi}{2}]^{2}$ and the
translation vector $\boldsymbol{\mathbf{v}}=(0,\frac{\pi}{2})^{\top}$. We
leave the recovering strategies of the whole region $[0,2\pi]^{2}$ and other
enhancement for future works.

######  Theorem 5.1.

For two different vectors
$\boldsymbol{\mathbf{x}}_{j}\in\mathbb{[}0,\frac{\pi}{2}]\times[\frac{\pi}{2},\pi],j=1,2$,
if $||\boldsymbol{\mathbf{x}}_{1}-\boldsymbol{\mathbf{x}}_{2}||_{1}\geq C$ for
a constant $C$, then

$\Big{|}{e^{i\boldsymbol{\mathbf{x}}_{1,1}}+e^{i\boldsymbol{\mathbf{x}}_{1,2}}-(e^{i\boldsymbol{\mathbf{x}}_{2,1}}+e^{i\boldsymbol{\mathbf{x}}_{2,2}})}\Big{|}\geq\frac{2C}{\pi^{2}}C.$

###### Proof.

We prove the lemma by considering the following two cases.
Case 1:
$0\leq\boldsymbol{\mathbf{x}}_{1,1}\leq\boldsymbol{\mathbf{x}}_{2,1}\leq\boldsymbol{\mathbf{x}}_{2,2}\leq\boldsymbol{\mathbf{x}}_{1,2}\leq\pi$.
In this case,

$\displaystyle\Big{|}{e^{i\boldsymbol{\mathbf{x}}_{1,1}}+e^{i\boldsymbol{\mathbf{x}}_{1,2}}-(e^{i\boldsymbol{\mathbf{x}}_{2,1}}+e^{i\boldsymbol{\mathbf{x}}_{2,2}})}\Big{|}\geq$
$\displaystyle\Big{|}{e^{i\boldsymbol{\mathbf{x}}_{2,1}}+e^{i\boldsymbol{\mathbf{x}}_{2,2}}}\Big{|}-\Big{|}{e^{i\boldsymbol{\mathbf{x}}_{1,1}}+e^{i\boldsymbol{\mathbf{x}}_{1,2}}}\Big{|}$
$\displaystyle\geq$ $\displaystyle
2\Big{(}\cos(\frac{\phi_{2}}{2})-\cos(\frac{\phi_{1}}{2})\Big{)},$

where
$\phi_{j}=\boldsymbol{\mathbf{x}}_{j,2}-\boldsymbol{\mathbf{x}}_{j,1},j=1,2$.
By the assumption of the theorem, we have $C\leq\phi_{1}-\phi_{2}\leq\pi$ and
$C\leq\phi_{1}+\phi_{2}\leq 2\pi$. Thus

$2\Big{(}\cos(\frac{\phi_{2}}{2})-\cos(\frac{\phi_{1}}{2})\Big{)}=4\sin(\frac{\phi_{1}+\phi_{2}}{4})\sin(\frac{\phi_{1}-\phi_{2}}{4})\geq
4\sin(\frac{C}{4})\sin(\frac{C}{4})\geq\frac{2C^{2}}{\pi^{2}}.$

where the last inequality uses
$\sin(\frac{C}{4})\geq\frac{2\sqrt{2}}{\pi}\frac{C}{4}$ for
$0\leq\frac{C}{4}\leq\frac{\pi}{4}$.
Case 2:
$0\leq\boldsymbol{\mathbf{x}}_{1,1}\leq\boldsymbol{\mathbf{x}}_{2,1}\leq\boldsymbol{\mathbf{x}}_{1,2}\leq\boldsymbol{\mathbf{x}}_{2,2}\leq\pi$.
Again, the idea is to calculate the angle between
$e^{i\boldsymbol{\mathbf{x}}_{1,1}}+e^{i\boldsymbol{\mathbf{x}}_{1,2}}$ and
$e^{i\boldsymbol{\mathbf{x}}_{2,1}}+e^{i\boldsymbol{\mathbf{x}}_{2,2}}$. By a
simple analysis of the angle relations between
$e^{i\boldsymbol{\mathbf{x}}_{1,1}},e^{i\boldsymbol{\mathbf{x}}_{1,2}},e^{i\boldsymbol{\mathbf{x}}_{2,1}},e^{i\boldsymbol{\mathbf{x}}_{2,2}},$
we obtain that the angle between
$e^{i\boldsymbol{\mathbf{x}}_{1,1}}+e^{i\boldsymbol{\mathbf{x}}_{1,2}}$ and
$e^{i\boldsymbol{\mathbf{x}}_{2,1}}+e^{i\boldsymbol{\mathbf{x}}_{2,2}}$ is
$\frac{\boldsymbol{\mathbf{x}}_{2,1}-\boldsymbol{\mathbf{x}}_{1,1}+\boldsymbol{\mathbf{x}}_{2,2}-\boldsymbol{\mathbf{x}}_{1,2}}{2}$
which is larger than $\frac{C}{2}$. Thus

$\Big{|}{e^{i\boldsymbol{\mathbf{x}}_{1,1}}+e^{i\boldsymbol{\mathbf{x}}_{1,2}}-(e^{i\boldsymbol{\mathbf{x}}_{2,1}}+e^{i\boldsymbol{\mathbf{x}}_{2,2}})}\Big{|}\geq\max\Big{(}\Big{|}{e^{i\boldsymbol{\mathbf{x}}_{1,1}}+e^{i\boldsymbol{\mathbf{x}}_{1,2}}}\Big{|},\Big{|}{e^{i\boldsymbol{\mathbf{x}}_{2,1}}+e^{i\boldsymbol{\mathbf{x}}_{2,2}}}\Big{|}\Big{)}\sin(\frac{C}{2}).$
(5.5)

We next claim that

$\max\Big{(}\Big{|}{e^{i\boldsymbol{\mathbf{x}}_{1,1}}+e^{i\boldsymbol{\mathbf{x}}_{1,2}}}\Big{|},\Big{|}{e^{i\boldsymbol{\mathbf{x}}_{2,1}}+e^{i\boldsymbol{\mathbf{x}}_{2,2}}}\Big{|}\Big{)}\geq
2\cos\Big{(}\frac{\pi-C/2}{2}\Big{)}.$

Otherwise,
$\boldsymbol{\mathbf{x}}_{1,2}-\boldsymbol{\mathbf{x}}_{1,1}>\pi-\frac{C}{2}$
and
$\boldsymbol{\mathbf{x}}_{2,2}-\boldsymbol{\mathbf{x}}_{2,1}>\pi-\frac{C}{2}$,
which is impossible when
$||\boldsymbol{\mathbf{x}}_{1}-\boldsymbol{\mathbf{x}}_{2}||_{1}\geq C$. Thus
the claim is proved. Together with (5.5), we arrive at

$\Big{|}{e^{i\boldsymbol{\mathbf{x}}_{1,1}}+e^{i\boldsymbol{\mathbf{x}}_{1,2}}-(e^{i\boldsymbol{\mathbf{x}}_{2,1}}+e^{i\boldsymbol{\mathbf{x}}_{2,2}})}\Big{|}\geq
2\sin(\frac{C}{4})\sin(\frac{C}{2})\geq\frac{2C^{2}}{\pi^{2}}.$

This completes the proof. ∎

Figure 5.1: Although the sources are well-separated, the direction components
of sources are closely spaced.

####  5.4.2 Phase transition and performance of Algorithm 4

Most of the conventional two-dimensional DOA algorithms consider multiple
snapshots of measurements from coherent or incoherent signals. Also, the noise
is usually assumed to be white Gaussian noise such that the expectation of the
covariance matrix of the measurement vector is a sum of two terms, where the
first term is from the correlation of the signals and the second one is the
noise correlation matrix. Based on this crucial observation, many algorithms
were derived to tackle the problem. Differently to the above model, we
consider recovering the source from a single measurement with deterministic
noise. Thus we do not compare the performance of our algorithm with those
algorithms with statistical model. We demonstrate the super-resolution
capacity of our algorithm for the single snapshot case by showing the phase
transition of the algorithm. We will derive a coordinate-combination-based
MUSIC algorithm for multiple snapshots case in a forthcoming work.

We now describe the numerical experiments for demonstrating the phase
transition phenomenon of our algorithm in terms of the SNR versus the super-
resolution factor. We fix $\Omega=10$ and consider three and four sources
separated by the minimum separation distance $D_{\min}$, i.e., $\min_{p\neq
q}||\boldsymbol{\mathbf{y}}_{p}-\boldsymbol{\mathbf{y}}_{q}||_{1}\geq
D_{\min}$. We perform 10000 random experiments (the randomness is in the
choice of $(D_{\min},\sigma,\boldsymbol{\mathbf{y}}_{j},a_{j})$ to recover the
sources using Algorithm 4. The reconstruction is viewed and recorded as
successful if the recovered source is in a $\frac{D_{\min}}{3}$-neighborhood
of the underlying source, otherwise it is unsuccessful; See Algorithm 5 for
the details of a single experiment. The results of the experiments are
summarized in Figure 4.2 which shows each successful and unsuccessfully
recovery with respective to the $\log(SRF)$ and $\log(SNR)$. It is observed
that there is a line with slope ($2n-1$) in the parameter space $\log(SRF)$
versus $\log(SNR)$ above which the source is stably reconstructed for every
realization. This phase transition phenomenon is exactly the one predicted by
our theoretical result in Theorems 2.2. It also manifests the efficiency of
Algorithm 4 as it can resolve the source in the regime where the source
separation distance is of the order of the computational resolution limit.

Input: Sources $\mu=\sum_{j=1}^{n}a_{j}\delta_{\boldsymbol{\mathbf{y}}_{j}}$;
Noise level $\sigma$;

Input: Measurements:
$\mathbf{Y}(\boldsymbol{\mathbf{\omega}}),\boldsymbol{\mathbf{\omega}}=[0,1,\cdots,\Omega]^{2}$;

1: $\text{Successnumber}=0$;

2: Input source number $n$ and measurement $\boldsymbol{\mathbf{Y}}$ to
Algorithm 4 and save the output as
$\boldsymbol{\mathbf{y}}_{1},\cdots,\boldsymbol{\mathbf{y}}_{n}$;

for  _each $1\leq j\leq n$_ do

Compute the error for the source location $\boldsymbol{\mathbf{y}}_{j}$:
$e_{j}:=\min_{\mathbf{\hat{y}}_{l},l=1,\cdots,n}||\mathbf{\hat{y}}_{l}-\boldsymbol{\mathbf{y}}_{j}||_{2}$;

The source location $\boldsymbol{\mathbf{y}}_{j}$ is recovered successfully if
$e_{j}<\frac{\min_{p\neq
j}||\boldsymbol{\mathbf{y}}_{p}-\boldsymbol{\mathbf{y}}_{j}||_{2}}{3};$ and
$\text{Successnumber}=\text{Successnumber}+1;$

if  _$\text{Successnumber}==n$_ then

Return Success

else

Return Fail

.

Algorithm 5 A single experiment

(a) Recovery success.

(b) Recovery success.

Figure 5.2: Plots of the successful and the unsuccessful location recoveries
by Algorithm 4 in terms of $\log(\frac{1}{\sigma})$ versus $\log(SRF)$. (a)
illustrates that locations of three point sources can be stably recovered if
$\log(\frac{1}{\sigma})$ is above a line of slope $5$ in the parameter space.
Conversely, for the same case, (b) shows that locations of four point sources
can be stably recovered if $\log(\frac{1}{\sigma})$ is above a line of slope
$7$ in the parameter space.

##  6 A nonlinear approximation theory in Vandermonde space

In this section, we introduce the main technique, a nonlinear approximation
theory in Vandermonde space[44, 43], that is used to deal with one-dimensional
super-resolution problems. In [44], we have derived the theory for real
numbers and in [43] for complex numbers on the unit circle. Here, we derive a
different theory for arbitrary bounded complex numbers, which are related to
the proofs of the main results of the paper.

For a given positive integer $s$ and $\omega\in\mathbb{C}$, we denote by

$\phi_{s}(\omega)=(1,\omega,\cdots,\omega^{s})^{\top}$ (6.1)

and call $\phi_{s}$ a Vandermonde vector. At the heart of the theory is the
following nonlinear approximation problem in the Vandermonde space

$\min_{\hat{a}_{j},\hat{d}_{j}\in\mathbb{R},|\hat{d}_{j}|\leq
d,j=1,\cdots,k}\Big{|}\Big{|}\sum_{j=1}^{k}\hat{a}_{j}\phi_{s}(\hat{d}_{j})-v\Big{|}\Big{|}_{2},$
(6.2)

where $v=\sum_{j=1}^{k+1}a_{j}\phi_{s}(d_{j})$ is a given vector. We shall
derive a sharp lower-bound for this problem. In addition, we shall also
investigate the stability of the approximation problem (6.2) for
$v=\sum_{j=1}^{k}a_{j}\phi_{s}(d_{j})$.

###  6.1 Notation and Preliminaries

We introduce some notation and preliminaries. We denote the Vandermonde matrix
by

$\displaystyle V_{s}(k)=\begin{pmatrix}1&\cdots&1\\\ d_{1}&\cdots&d_{k}\\\
\vdots&\ddots&\vdots\\\
d_{1}^{s}&\cdots&d_{k}^{s}\end{pmatrix}=\Big{(}\phi_{s}(d_{1})\ \
\phi_{s}(d_{2})\ \ \cdots\ \ \phi_{s}(d_{k})\Big{)}.$ (6.3)

For a real matrix or a vector $A$, we denote by $A^{\top}$ its transpose and
by $A^{*}$ its conjugate transpose.

We first present some basic properties of Vandermonde matrices.

######  Lemma 6.1.

For $k$ distinct complex numbers $d_{j}$’s, we have

$||V_{k-1}(k)^{-1}||_{\infty}\leq\max_{1\leq i\leq k}\Pi_{1\leq p\leq k,p\neq
i}\frac{1+|d_{p}|}{|d_{i}-d_{p}|},$

where $V_{k-1}(k)$ is the Vandermonde matrix $V_{k-1}(k)$ defined as in (6.3).

###### Proof.

See Theorem 1 in [26].∎

As a consequence, we directly have the following corollary.

######  Corollary 6.2.

Let $d_{\min}=\min_{i\neq j}|d_{i}-d_{j}|$ and assume that
$\max_{i=1,\cdots,k}|d_{i}|\leq d$. Then

$||V_{k-1}(k)^{-1}||_{\infty}\leq\frac{(1+d)^{k-1}}{(d_{\min})^{k-1}}.$

######  Lemma 6.3.

For distinct $d_{1},\cdots,d_{k}\in\mathbb{C}$, define the Vandermonde
matrices $V_{k-1}(k),V_{s}(k)$ as in (6.3) with $s\geq k-1$. Then the
following estimate on their singular values holds:

$\frac{1}{\sqrt{k}}\frac{1}{||V_{k-1}(k)^{-1}||_{\infty}}\leq\frac{1}{||V_{k-1}(k)^{-1}||_{2}}\leq\sigma_{\min}(V_{k-1}(k))\leq\sigma_{\min}(V_{s}(k)).$

###### Proof.

The result holds by using properties of matrix norms. ∎

Denote by

$S_{1k}^{j}:=\Big{\\{}\\{\tau_{1},\cdots,\tau_{j}\\}:\text{$\tau_{p}\in\\{1,\cdots,k\\},p=1,\cdots,j$
and $\tau_{p}\neq\tau_{q}$, for $p\neq q$}\Big{\\}}.$

Note that there is no order in $\\{\tau_{1},\cdots,\tau_{j}\\}$, i.e.,
$\\{1,2\\}$ and $\\{2,1\\}$ are the same sets. We then have the following
decomposition of the Vandermonde matrix.

######  Proposition 6.4.

The Vandermonde matrix $V_{k}(k)$ defined as in (6.3) can be reduced to the
following form by using elementary column-addition operations, i.e.,

$\displaystyle V_{k}(k)G(1)\cdots G(k-1)DQ(1)\cdots
Q(k-1)=\begin{pmatrix}1&0&\cdots&0\\\ 0&1&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots\\\ 0&0&\cdots&1\\\
v_{(k+1)1}&v_{(k+1)2}&\cdots&v_{(k+1)k}\end{pmatrix},$ (6.4)

where $G(1),\cdots,G(k-1),Q(1),\cdots,Q(k-1)$ are elementary column-addition
matrices,

$D=\text{diag}(1,\frac{1}{(d_{2}-d_{1})},\cdots,\frac{1}{\Pi_{p=1}^{k-1}(d_{k}-d_{p})})$

and

$v_{(k+1)j}=(-1)^{k-j}\sum_{\\{\tau_{1},\cdots,\tau_{k+1-j}\\}\in
S_{1k}^{k+1-j}}d_{\tau_{1}}\cdots d_{\tau_{k+1-j}}.$ (6.5)

###### Proof.

See Appendix B in [44]. ∎

######  Lemma 6.5.

For an $s\times k$ complex matrix $A$ of rank $k$ with $s>k$, let $V$ be the
space spanned by columns of $A$ and $V^{\perp}$ be the orthogonal complement
of $V$. Denote by $P_{V^{\perp}}$ the orthogonal projection to $V^{\perp}$,
and set $D=(A,v)$. We have

$\min_{a\in\mathbb{C}^{k}}||Aa-v||_{2}=||P_{V^{\perp}}(v)||_{2}=\sqrt{\frac{\det(D^{*}D)}{\det(A^{*}A)}}.$

###### Proof.

See Lemma 1 in [43]. ∎

######  Lemma 6.6.

We have

$\sqrt{\frac{\det(V_{k}(k)^{*}V_{k}(k))}{\det(V_{k-1}(k)^{*}V_{k-1}(k))}}=\sqrt{\sum_{j=0}^{k}|v_{j}|^{2}},$
(6.6)

where $V_{s}(k)$ is defined as in (6.3) and
$v_{j}=\sum_{\\{\tau_{1},\cdots,\tau_{j}\\}\in S_{1k}^{j}}d_{\tau_{1}}\cdots
d_{\tau_{j}}$. Especially, if $|d_{j}|<d,j=1,\cdots,k$, then

$\sqrt{\frac{\det(V_{k}(k)^{*}V_{k}(k))}{\det(V_{k-1}(k)^{*}V_{k-1}(k))}}\leq(1+d)^{k}.$
(6.7)

###### Proof.

Note that in Proposition 6.4, all the elementary column-addition matrices have
unit determinant. As a result,
$\det(V_{k}(k)^{*}V_{k}(k))=\frac{\det(F^{*}F)}{\det(D^{*}D)}$, where $F$ is
the matrix in the right-hand side of (6.4), and $D$ is the diagonal matrix in
Proposition 6.4. A direct calculation shows that
$\det(F^{*}F)=\sum_{j=0}^{k}|v_{j}|^{2}$, where we use (6.5). On the other
hand, $V_{k-1}(k)$ is a standard Vandermonde matrix and we have
$\det(V_{k-1}(k)^{*}V_{k-1}(k))=\frac{1}{\det(D^{*}D)}$. Combining these
results, (6.6) follows. The last statement can be derived from (6.6) and the
estimate that

$\sqrt{\sum_{j=0}^{k}|v_{j}|^{2}}\leq\sum_{j=0}^{k}|v_{j}|\leq\sum_{j=0}^{k}\begin{pmatrix}k\\\
j\end{pmatrix}d^{j}=(1+d)^{k}.$

∎

For reader’s convenience, we finally present two auxiliary lemmas. For
positive integers $p,q$ and complex numbers
$z_{1},\cdots,z_{p},\hat{z}_{1},\cdots,\hat{z}_{q}$, we define

$\eta_{p,q}(z_{1},\cdots,z_{p},\hat{z}_{1},\cdots,\hat{z}_{q})=\begin{pmatrix}|z_{1}-\hat{z}_{1}|\cdots|z_{1}-\hat{z}_{q}|\\\
|z_{2}-\hat{z}_{1}|\cdots|z_{2}-\hat{z}_{q}|\\\ \vdots\\\
|z_{p}-\hat{z}_{1}|\cdots|z_{p}-\hat{z}_{q}|\end{pmatrix}.$ (6.8)

The following two properties of $\eta_{p,q}$ hold.

######  Lemma 6.7.

For complex numbers $d_{j},\hat{d}_{j}$’s, we have the following estimate

$\Big{|}\Big{|}\eta_{k+1,k}(d_{1},\cdots,d_{k+1},\hat{d}_{1},\cdots,\hat{d}_{k})\Big{|}\Big{|}_{\infty}\geq(\frac{d_{\min}}{2})^{k},$

where $d_{\min}=\min_{j\neq p}|d_{j}-d_{p}|$ and
$\eta_{k+1,k}(d_{1},\cdots,d_{k+1},\hat{d}_{1},\cdots,\hat{d}_{k})$ is defined
as in (6.8).

###### Proof.

Because we have $k+1$ $d_{j}$’s and only $k$ $\hat{d}_{j}$’s, there must exist
one $d_{j_{0}}$ so that

$|d_{j_{0}}-\hat{d}_{j}|\geq\frac{d_{\min}}{2},\quad j=1,\cdots,k.$

Then the estimate in the lemma follows. ∎

######  Lemma 6.8.

Let $d_{j},\hat{d}_{j}\in\mathbb{C},j=1,\cdots,k$ satisfy
$|d_{j}|,|\hat{d}_{j}|\leq d$. Assume that

$||\eta_{k,k}(d_{1},\cdots,d_{k},\hat{d}_{1},\cdots,\hat{d}_{k})||_{\infty}<\epsilon,$
(6.9)

where $\eta_{k,k}(\cdots)$ is defined as in (6.8), and that

$d_{\min}=\min_{p\neq q}|d_{p}-d_{j}|\geq 2\epsilon^{\frac{1}{k}}.$ (6.10)

Then after reordering $d_{j}$’s, we have

$\Big{|}{\hat{d}_{j}-d_{j}}\Big{|}<\frac{d_{\min}}{2},\quad j=1,\cdots,k,$
(6.11)

and moreover

$\Big{|}{\hat{d}_{j}-d_{j}}\Big{|}\leq\Big{(}\frac{2}{d_{\min}}\Big{)}^{k-1}\epsilon,\quad
j=1,\cdots,k.$ (6.12)

###### Proof.

See Appendix A. ∎

###  6.2 Lower-bound for the approximation problem (6.2)

In this section, we derive a lower-bound for the nonlinear approximation
problem (6.2). We first consider the special case when $v$ is a Vandermonde
vector.

######  Theorem 6.1.

Let $k\geq 1$ and $\hat{d}_{1},\cdots,\hat{d}_{k}$ be $k$ distinct complex
numbers with $|\hat{d}_{j}|\leq\hat{d},1\leq j\leq k$. Define
$A:=\big{(}\phi_{k}(\hat{d}_{1}),\cdots,\phi_{k}(\hat{d}_{k})\big{)}$, where
$\phi_{k}(\hat{d}_{j})$’s are defined as in (6.1). Let $V$ be the
$k$-dimensional space spanned by the column vectors of $A$, and let
$V^{\perp}$ be the one-dimensional orthogonal complement of $V$ in
$\mathbb{C}^{k+1}$. Let $P_{V^{\perp}}$ be the orthogonal projection onto
$V^{\perp}$ in $\mathbb{C}^{k+1}$. Then we have

$\min_{a\in\mathbb{C}^{k}}\Big{|}\Big{|}{Aa-\phi_{k}(x)}\Big{|}\Big{|}_{2}=\Big{|}\Big{|}{P_{V^{\perp}}(\phi_{k}(x))}\Big{|}\Big{|}_{2}=\Big{|}{v^{*}\phi_{k}(x)}\Big{|}\geq\frac{1}{(1+\hat{d})^{k}}\Big{|}{\Pi_{j=1}^{k}(x-\hat{d}_{j})}\Big{|},$

where $v$ is a unit vector in $V^{\perp}$ and $v^{*}$ is its conjugate
transpose.

###### Proof.

By Lemma 6.5, it follows that

$\min_{a\in\mathbb{C}^{k}}\Big{|}\Big{|}{Aa-\phi_{k}(x)}\Big{|}\Big{|}_{2}=\sqrt{\frac{\text{det}(D^{*}D)}{\text{det}(A^{*}A)}},$

where
$D=\big{(}\phi_{k}(\hat{d}_{1}),\cdots,\phi_{k}(\hat{d}_{k}),\phi_{k}(x)\big{)}$.
Denote
$\tilde{A}=\big{(}\phi_{k-1}(\hat{d}_{1}),\cdots,\phi_{k-1}(\hat{d}_{k})\big{)}$.
By (6.7), we have

$\sqrt{\frac{\text{det}(A^{*}A)}{\text{det}(\tilde{A}^{*}\tilde{A})}}\leq(1+\hat{d})^{k}.$

Therefore,

$\min_{a\in\mathbb{C}^{k}}\Big{|}\Big{|}{Aa-\phi_{k}(x)}\Big{|}\Big{|}_{2}\geq\frac{1}{(1+\hat{d})^{k}}\sqrt{\frac{\text{det}(D^{*}D)}{\text{det}(\tilde{A}^{*}\tilde{A})}}.$

Note that $D$ and $\tilde{A}$ are square Vandermonde matrices. We can use the
determinant formula to derive that

$\min_{a\in\mathbb{C}^{k}}\Big{|}\Big{|}{Aa-\phi_{k}(x)}\Big{|}\Big{|}_{2}\geq\frac{1}{(1+\hat{d})^{k}}\frac{|\Pi_{1\leq
t<p\leq k}(\hat{d}_{t}-\hat{d}_{p})\Pi_{q=1}^{k}(x-\hat{d}_{q})|}{|\Pi_{1\leq
t<p\leq
k}(\hat{d}_{t}-\hat{d}_{p})|}=\frac{1}{(1+\hat{d})^{k}}|\Pi_{j=1}^{k}(x-\hat{d}_{j})|.$

This completes the proof of the theorem. ∎

We now consider the approximation problem (6.2) for the general case when $v$
is a linear combination of Vandermonde vectors.

######  Theorem 6.2.

Let $k\geq 1$. Assume $(k+1)$ different complex numbers
$d_{j}\in\mathbb{C},j=1,\cdots,k+1$ with $|d_{j}|\leq d$ and $(k+1)$
$a_{j}\in\mathbb{C}$ with $|a_{j}|\geq m_{\min}$. Let $d_{\min}:=\min_{j\neq
p}|d_{j}-d_{p}|$. For $q\leq k$, let
$\hat{a}(q)=(\hat{a}_{1},\hat{a}_{2},\cdots,\hat{a}_{q})^{\top},a=(a_{1},a_{2},\cdots,a_{k+1})^{\top}$,
and

$\hat{A}(q)=\big{(}\phi_{2k}(\hat{d}_{1}),\ \cdots,\
\phi_{2k}(\hat{d}_{q})\big{)},\ A=\big{(}\phi_{2k}(d_{1}),\ \cdots,\
\phi_{2k}(d_{k+1})\big{)},$

where $\phi_{2k}(z)$ is defined as in (6.1). Then

$\displaystyle\min_{\hat{a}_{p},\hat{d}_{p}\in\mathbb{C},|\hat{d}_{p}|\leq\hat{d},p=1,\cdots,q}||\hat{A}(q)\hat{a}(q)-Aa||_{2}\geq\frac{m_{\min}(d_{\min})^{2k}}{2^{k}(1+d)^{k}(1+\hat{d})^{k}}.$

###### Proof.

Step 1. Note that for $q<k$, we have

$\min_{\hat{a}_{p},\hat{d}_{p}\in\mathbb{C},|\hat{d}_{p}|\leq
d,p=1,\cdots,q}||\hat{A}(q)\hat{a}(q)-Aa||_{2}\geq\min_{\hat{a}_{p},\hat{d}_{p}\in\mathbb{C},|\hat{d}_{p}|\leq
d,p=1,\cdots,k}||\hat{A}(k)\hat{a}(k)-Aa||_{2}.$

Hence we need only to consider the case when $q=k$. It then suffices to show
that for any given
$\hat{d}_{j}\in\mathbb{C},|\hat{d}_{j}|\leq\hat{d},j=1,\cdots,k$, the
following holds

$\min_{\hat{a}_{p}\in\mathbb{C},p=1,\cdots,k}||\hat{A}(k)\hat{a}(k)-Aa||_{2}\geq\frac{m_{\min}(d_{\min})^{2k}}{2^{k}(1+d)^{k}(1+\hat{d})^{k}}.$
(6.13)

So we fix $\hat{d}_{1},\cdots,\hat{d}_{k}$ in our subsequent argument.
Step 2. For $l=0,\cdots,k$, we define the following partial matrices

$\hat{A}_{l}=\left(\begin{array}[]{ccc}\hat{d}_{1}^{l}&\cdots&\hat{d}_{k}^{l}\\\
\hat{d}_{1}^{l+1}&\cdots&\hat{d}_{k}^{l+1}\\\ \vdots&\vdots&\vdots\\\
\hat{d}_{1}^{l+k}&\cdots&\hat{d}_{k}^{l+k}\end{array}\right),\quad
A_{l}=\left(\begin{array}[]{ccc}(d_{1})^{l}&\cdots&(d_{k+1})^{l}\\\
(d_{1})^{l+1}&\cdots&(d_{k+1})^{l+1}\\\ \vdots&\vdots&\vdots\\\
(d_{1})^{l+k}&\cdots&(d_{k+1})^{l+k}\end{array}\right).$

It is clear that for all $l$,

$\min_{\hat{a}(k)\in\mathbb{C}^{k}}||\hat{A}(k)\hat{a}(k)-Aa||_{2}\geq\min_{\hat{a}\in\mathbb{C}^{k}}||\hat{A}_{l}\hat{a}-A_{l}a||_{2}.$
(6.14)

Step 3. For each $l$, observe that
$\hat{A}_{l}=\hat{A}_{0}\text{diag}(\hat{d}_{1}^{l},\ \cdots,\
\hat{d}_{k}^{l}),A_{l}=A_{0}\text{diag}(d_{1}^{l},\cdots,d_{k+1}^{l})$, and
thus

$\min_{\hat{a}\in\mathbb{C}^{k}}||\hat{A}_{l}\hat{a}-A_{l}a||_{2}\geq\min_{\hat{\alpha}_{l}\in\mathbb{C}^{k}}||A_{0}\hat{\alpha}_{l}-A_{0}\alpha_{l}||_{2},$
(6.15)

where
$\alpha_{l}=\left(a_{1}(d_{1})^{l},\cdots,a_{k+1}(d_{k+1})^{l}\right)^{\top}$.
Let $V$ be the space spanned by the column vectors of $A_{0}$. Then the
dimension of $V$ is $k$, and the dimension of $V^{\perp}$, the orthogonal
complement of $V$ in $\mathbb{C}^{k+1}$, is one. Let $P_{V^{\perp}}$ be the
orthogonal projection onto $V^{\perp}$. Note that
$||P_{V^{\perp}}u||_{2}=|v^{*}u|$ for $u\in\mathbb{R}^{k+1}$, where $v$ is a
unit vector in $V^{\perp}$ and $v^{*}$ is its conjugate transpose. We have

$\displaystyle\min_{\hat{\alpha}_{l}\in\mathbb{C}^{k}}||\hat{A}_{0}\hat{\alpha}_{l}-A_{0}\alpha_{l}||_{2}=||P_{V^{\perp}}(A_{0}\alpha_{l})||_{2}=|v^{*}A_{0}\alpha_{l}|=\Big{|}\sum_{j=1}^{k+1}a_{j}(d_{j})^{l}v^{*}\phi_{k}(d_{j})\Big{|}=|\beta_{l}|,$
(6.16)

where

$\beta_{l}=\sum_{j=1}^{k+1}a_{j}(d_{j})^{l}v^{*}\phi_{k}(d_{j}),\quad\text{for
$l=0,1,\cdots,k.$}$

Step 4. Denote $\beta=(\beta_{0},\cdots,\beta_{k})^{\top}$. We have
$B\hat{\eta}=\beta$, where

$B=\left(\begin{array}[]{cccc}a_{1}&a_{2}&\cdots&a_{k+1}\\\
a_{1}d_{1}&a_{2}d_{2}&\cdots&a_{k+1}d_{k+1}\\\ \vdots&\vdots&\vdots&\vdots\\\
a_{1}(d_{1})^{k}&a_{2}(d_{2})^{k}&\cdots&a_{k+1}(d_{k+1})^{k}\end{array}\right),\quad\hat{\eta}=\left(\begin{array}[]{c}v^{*}\phi_{k}(d_{1})\\\
v^{*}\phi_{k}(d_{2})\\\ \vdots\\\ v^{*}\phi_{k}(d_{k+1})\end{array}\right).$

Corollary 6.2 yields

$\displaystyle||\hat{\eta}||_{\infty}=||B^{-1}\beta||_{\infty}\leq||B^{-1}||_{\infty}||\beta||_{\infty}\leq\frac{(1+d)^{k}}{m_{\min}(d_{\min})^{k}}||\beta||_{\infty}.$

On the other hand, applying Theorem 6.1 to each term $|v^{*}\phi_{k}(d_{j})|$,
$j=1,2,\cdots k+1$, we have

$||\hat{\eta}||_{\infty}\geq\frac{1}{(1+\hat{d})^{k}}||\eta_{k+1,k}(d_{1},\cdots,d_{k+1},\hat{d}_{1},\cdots,\hat{d}_{k})||_{\infty},$

where $\eta_{k+1,k}(\cdots)$ is defined as in (6.8). Combining this inequality
with Lemma 6.7, we get

$||\hat{\eta}||_{\infty}\geq\frac{(d_{\min})^{k}}{2^{k}(1+\hat{d})^{k}}.$

Then it follows that

$||\beta||_{\infty}\geq\frac{m_{\min}(d_{\min})^{2k}}{2^{k}(1+d)^{k}(1+\hat{d})^{k}}.$

Therefore, recalling (6.14)–(6.16), we arrive at

$\min_{\hat{a}(k)\in\mathbb{C}^{k}}||\hat{A}(k)\hat{a}(k)-Aa||_{2}\geq\max_{0\leq
l\leq
k}\min_{\hat{a}\in\mathbb{C}^{k}}||\hat{A}_{l}\hat{a}-A_{l}a||_{2}=\max_{0\leq
l\leq
k}|\beta_{l}|=||\beta||_{\infty}\geq\frac{m_{\min}(d_{\min})^{2k}}{2^{k}(1+d)^{k}(1+\hat{d})^{k}}.$

This proves (6.13) and hence the theorem. ∎

###  6.3 Stability of the approximation problem (6.2)

In the section we present a stability result for the approximation problem
(6.2).

######  Theorem 6.3.

Let $k\geq 1$. Assume $k$ different complex numbers
$d_{j}\in\mathbb{C},j=1,\cdots,k$ with $|d_{j}|\leq d$ and $k$
$a_{j}\in\mathbb{C}$ with $|a_{j}|\geq m_{\min}$. Let $d_{\min}:=\min_{p\neq
q}|d_{p}-d_{q}|$. Assume that $\hat{d}_{j}\in\mathbb{C},j=1,\cdots,k$ with
$|\hat{d}_{j}|\leq d$ satisfy

$||\hat{A}\hat{a}-Aa||_{2}<\sigma,$

where $\hat{a}=(\hat{a}_{1},\cdots,\hat{a}_{k})^{\top}$,
$a=(a_{1},\cdots,a_{k})^{\top}$, and

$\hat{A}=\big{(}\phi_{2k-1}(\hat{d}_{1}),\ \cdots,\
\phi_{2k-1}(\hat{d}_{k})\big{)},\ A=\big{(}\phi_{2k-1}(d_{1}),\ \cdots,\
\phi_{2k-1}(d_{k})\big{)}.$

Then

$\Big{|}\Big{|}\eta_{k,k}(d_{1},\cdots,d_{k},\hat{d}_{1},\cdots,\hat{d}_{k})\Big{|}\Big{|}_{\infty}<\frac{(1+d)^{2k-1}}{d_{\min}^{k-1}}\frac{\sigma}{m_{\min}}.$

###### Proof.

Since $||\hat{A}\hat{a}-Aa||_{2}<\sigma$, we have

$\min_{\hat{\alpha}\in\mathbb{C}^{k}}||\hat{A}\hat{\alpha}-Aa||_{2}<\sigma,$

and hence

$\max_{0\leq l\leq
k-1}\min_{\hat{\alpha}\in\mathbb{C}^{k}}||\hat{A}_{l}\hat{\alpha}-A_{l}a||_{2}\leq\min_{\hat{\alpha}\in\mathbb{C}^{k}}||\hat{A}\hat{\alpha}-Aa||_{2}<\sigma,$
(6.17)

where

$\hat{A}_{l}=\left(\begin{array}[]{ccc}\hat{d}_{1}^{l}&\cdots&\hat{d}_{k}^{l}\\\
\hat{d}_{1}^{l+1}&\cdots&\hat{d}_{k}^{l+1}\\\ \vdots&\vdots&\vdots\\\
\hat{d}_{1}^{l+k}&\cdots&\hat{d}_{k}^{l+k}\end{array}\right),\quad
A_{l}=\left(\begin{array}[]{ccc}d_{1}^{l}&\cdots&d_{k}^{l}\\\
d_{1}^{l+1}&\cdots&d_{k}^{l+1}\\\ \vdots&\vdots&\vdots\\\
d_{1}^{l+k}&\cdots&d_{k}^{l+k}\end{array}\right).$

For each $l$, from the decomposition
$\hat{A}_{l}=\hat{A}_{0}\text{diag}(\hat{d}_{1}^{l},\ \cdots,\
\hat{d}_{k}^{l}),A_{l}=A_{0}\text{diag}((d_{1})^{l},\cdots,(d_{k})^{l})$, we
get

$\min_{\hat{\alpha}\in\mathbb{C}^{k}}||\hat{A}_{l}\hat{\alpha}-A_{l}a||_{2}\geq\min_{\hat{\alpha}_{l}\in\mathbb{C}^{k}}||\hat{A}_{0}\hat{\alpha}_{l}-A_{0}\alpha_{l}||_{2},$
(6.18)

where $\alpha_{l}=(a_{1}(d_{1})^{l},\cdots,a_{k}(d_{k})^{l})^{\top}$. Let $V$
be the space spanned by the column vectors of $\hat{A}_{0}$. Then the
dimension of $V$ is $k$, and $V^{\perp}$, the orthogonal complement of $V$ in
$\mathbb{C}^{k+1}$ is of dimension one. We let $v$ be a unit vector in
$V^{\perp}$ and let $P_{V^{\perp}}$ be the orthogonal projection onto
$V^{\perp}$. Similarly to (6.16), we have

$\displaystyle\min_{\hat{\alpha}_{l}\in\mathbb{C}^{k}}||\hat{A}_{0}\hat{\alpha}_{l}-A_{0}\alpha_{l}||_{2}=||P_{V^{\perp}}(A_{0}\alpha_{l})||_{2}=|v^{*}A_{0}\alpha_{l}|=\Big{|}\sum_{j=1}^{k}a_{j}(d_{j})^{l}v^{*}\phi_{k}(d_{j})\Big{|}=|\beta_{l}|,$
(6.19)

where $\beta_{l}=\sum_{j=1}^{k}a_{j}(d_{j})^{l}v^{*}\phi_{k}(d_{j})$. Let
$\beta=(\beta_{0},\cdots,\beta_{k-1})^{\top}$. Moreover, similarly to Step 4
in the proof of Theorem 6.2, we have

$\Big{|}\Big{|}\eta_{k,k}(d_{1},\cdots,d_{k},\hat{d}_{1},\cdots,\hat{d}_{k})\Big{|}\Big{|}_{\infty}\leq\frac{(1+d)^{2k-1}}{m_{\min}(d_{\min})^{k-1}}||\beta||_{\infty}.$

On the other hand, (6.17)–(6.19) indicate that $||\beta||_{\infty}<\sigma$.
Hence, we obtain that

$\Big{|}\Big{|}\eta_{k,k}(d_{1},\cdots,d_{k},\hat{d}_{1},\cdots,\hat{d}_{k})\Big{|}\Big{|}_{\infty}\leq\frac{(1+d)^{2k-1}}{(d_{\min})^{k-1}}\frac{\sigma}{m_{\min}}.$

This completes the proof. ∎

##  7 Conclusions and future works

In this paper, we have improved the estimates of resolution limits in two-
dimensional super-resolution problems. We also theoretically demonstrate the
optimal performance of a sparsity-promoting algorithm. Leveraging the new
techniques in the proof, we have proposed a coordinate-combination-based model
order detection algorithm and a coordinate-combination-based MUSIC algorithm
for DOA estimation in two dimensions. The superiority of the introduced
algorithms were demonstrated both theoretically or numerically.

Our work is also a start of many new topics. Firstly, one could extend the
techniques to three- and $k$-dimensional spaces to improve the resolution
estimates in higher dimensional super-resolution problems. Secondly, the idea
of coordinate-combination could inspire new algorithms for two-dimensional DOA
estimations in the case of multiple snapshots. These works will be presented
in a near future.

##  Appendix A Proof of Lemma 6.8

###### Proof.

Step 1. We claim that for each $\hat{d}_{p},1\leq p\leq k$, there exists one
$d_{j}$ such that $|\hat{d}_{p}-d_{j}|<\frac{d_{\min}}{2}$. By contradiction,
suppose that there exists $p_{0}$ such that
$|d_{j}-\hat{d}_{p_{0}}|\geq\frac{d_{\min}}{2}$ for all $1\leq j\leq k$.
Observe that

$\displaystyle\eta_{k,k}(d_{1},\cdots,d_{k},\hat{d}_{1},\cdots,\hat{d}_{k})$
$\displaystyle=$
$\displaystyle\text{diag}\left(|d_{1}-\hat{d}_{p_{0}}|,\cdots,|d_{k}-\hat{d}_{p_{0}}|\right)\eta_{k,k-1}(d_{1},\cdots,d_{k},\hat{d}_{1},\cdots,\hat{d}_{p_{0}-1},\hat{d}_{p_{0}+1},\cdots,\hat{d}_{k}).$

We write

$\eta_{k,k}=\eta_{k,k}(d_{1},\cdots,d_{k},\hat{d}_{1},\cdots,\hat{d}_{k})\quad\mbox{and}\quad\eta_{k,k-1}=\eta_{k,k-1}(d_{1},\cdots,d_{k},\hat{d}_{1},\cdots,\hat{d}_{p_{0}-1},\hat{d}_{p_{0}+1},\cdots,\hat{d}_{k}).$

Using Lemma 6.7, we have

$||\eta_{k,k}||_{\infty}\geq\frac{d_{\min}}{2}||\eta_{k,k-1}||_{\infty}\geq\Big{(}\frac{d_{\min}}{2}\Big{)}^{k}\geq\epsilon,$

where we have used (6.10) in the last inequality above. This contradicts (6.9)
and hence proves our claim.
Step 2. We claim that for each $d_{j},1\leq j\leq k$, there exists one and
only one $\hat{d}_{p}$ such that

$|d_{j}-\hat{d}_{p}|<\frac{d_{\min}}{2}.$

It suffices to show that for each $d_{j},1\leq j\leq k$, there is only one
$\hat{d}_{p}$ such that $|d_{j}-\hat{d}_{p}|<\frac{d_{\min}}{2}$. By
contradiction, suppose that there exist $p_{1},p_{2},$ and $j_{0}$ such that
$|d_{j_{0}}-\hat{d}_{p_{1}}|<\frac{d_{\min}}{2},|d_{j_{0}}-\hat{d}_{p_{2}}|<\frac{d_{\min}}{2}$.
Then for all $j\neq j_{0}$, we have

$\Big{|}(d_{j}-\hat{d}_{p_{1}})(d_{j}-\hat{d}_{p_{2}})\Big{|}\geq\frac{(d_{\min})^{2}}{4}.$
(A.1)

Similarly to the argument in Step 1, we separate the factors involving
$\hat{d}_{p_{1}},\hat{d}_{p_{2}},d_{j_{0}}$ from $\eta_{k,k}$ and consider

$\eta_{k-1,k-2}=\eta_{k-1,k-2}(d_{1},\cdots,d_{j_{0}-1},d_{j_{0}+1},\cdots,d_{k},\hat{d}_{1},\cdots,\hat{d}_{p_{1}-1},\hat{d}_{p_{1}+1},\cdots,\hat{d}_{p_{2}-1},\hat{d}_{p_{2}+1},\cdots,\hat{d}_{k}).$

Note that the components of $\eta_{k-1,k-2}$ differ from those of $\eta_{k,k}$
only by the factors $|(d_{j}-\hat{d}_{p_{1}})(d_{j}-\hat{d}_{p_{2}})|$ for
$j=1,\cdots,j_{0}-1,j_{0}+1,\cdots,k$. We can show that

$||\eta_{k,k}||_{\infty}\geq\frac{(d_{\min})^{2}}{4}||\eta_{k-1,k-2}||_{\infty}\geq\epsilon,$

where we have used Lemma 6.7 and (6.10) for establishing the last inequality
above. This contradicts (6.9) and hence proves our claim.
Step 3. By the result in Step 2, we can reorder $\hat{d}_{j}$’s to get

$|\hat{d}_{j}-d_{j}|<\frac{d_{\min}}{2},\quad j=1,\cdots,k.$

We now prove (6.12). It is clear that
$|\hat{d}_{p}-d_{j}|>\frac{d_{\min}}{2},p\neq j$. Thus

$|(d_{j}-\hat{d}_{1})\cdots(d_{j}-\hat{d}_{k})|>|d_{j}-d_{j}|(\frac{d_{\min}}{2})^{k-1},\quad
j=1,2,\cdots,k.$ (A.2)

Further, we get

$|d_{j}-\hat{d}_{j}|<\Big{(}\frac{2}{d_{\min}}\Big{)}^{k-1}||\eta_{k,k}||_{\infty}\leq\Big{(}\frac{2}{d_{\min}}\Big{)}^{k-1}\epsilon,\quad
j=1,2,\cdots,k.$

This completes the proof of the lemma. ∎

## References

  * [1] Ernst Abbe.  Beiträge zur theorie des mikroskops und der mikroskopischen wahrnehmung.  Archiv für mikroskopische Anatomie, 9(1):413–468, 1873.
  * [2] Hirotogu Akaike.  Information theory and an extension of the maximum likelihood principle.  In Selected papers of hirotugu akaike, pages 199–213. Springer, 1998.
  * [3] Hirotugu Akaike.  A new look at the statistical model identification.  In Selected Papers of Hirotugu Akaike, pages 215–222. Springer, 1974.
  * [4] Andrey Akinshin, Dmitry Batenkov, and Yosef Yomdin.  Accuracy of spike-train fourier reconstruction for colliding nodes.  In 2015 International Conference on Sampling Theory and Applications (SampTA), pages 617–621. IEEE, 2015.
  * [5] Habib Ammari.  An introduction to mathematics of emerging biomedical imaging.  Mathematics $\&$ Applications.
  * [6] Habib Ammari, Josselin Garnier, Wenjia Jing, Hyeonbae Kang, Mikyoung Lim, Knut Sølna, and Han Wang.  Mathematical and statistical methods for multistatic imaging.  Lecture Notes in Mathematics.
  * [7] Jean-Marc Azais, Yohann De Castro, and Fabrice Gamboa.  Spike detection from inaccurate samplings.  Applied and Computational Harmonic Analysis, 38(2):177–195, 2015\.
  * [8] Dmitry Batenkov, Laurent Demanet, Gil Goldman, and Yosef Yomdin.  Conditioning of partial nonuniform fourier matrices with clustered nodes.  SIAM Journal on Matrix Analysis and Applications, 41(1):199–220, 2020.
  * [9] Dmitry Batenkov, Gil Goldman, and Yosef Yomdin.  Super-resolution of near-colliding point sources.  Information and Inference: A Journal of the IMA, 05 2020.  iaaa005.
  * [10] E. Betzig, G.H. Patterson, R. Sougrat, O.W. Lindwasser, S. Olenych, J.S. Bonifacino, M.W. Davidson, J. Lippincott-Schwartz, and H.F. Hess.  Imaging intracellular fluorescent proteins at nanometer resolution.  Science, 313:1642–1645, 2006.
  * [11] Emmanuel J. Candès and Carlos Fernandez-Granda.  Towards a mathematical theory of super-resolution.  Commun. Pure Appl. Math., 67(6):906–956, 2014.
  * [12] Sitan Chen and Ankur Moitra.  Algorithmic foundations for the diffraction limit.  In Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing, pages 490–503, 2021.
  * [13] Weiguo Chen, Kon Max Wong, and James P Reilly.  Detection of the number of signals: A predicted eigen-threshold approach.  IEEE Transactions on Signal Processing, 39(5):1088–1098, 1991.
  * [14] Y-H Chen and C-H Chen.  Direction-of-arrival and frequency estimations for narrowband sources using two single rotation invariance algorithms with the marked subspace.  In IEE Proceedings F (Radar and Signal Processing), volume 139, pages 297–300. IET, 1992.
  * [15] Edward AK Cohen, Anish V Abraham, Sreevidhya Ramakrishnan, and Raimund J Ober.  Resolution limit of image analysis algorithms.  Nature Communications, 10(1):1–12, 2019.
  * [16] Maxime Ferreira Da Costa and Yuejie Chi.  On the stable resolution limit of total variation regularization for spike deconvolution.  IEEE Transactions on Information Theory, 66(11):7237–7252, 2020\.
  * [17] Geoffrey De Villiers and E Roy Pike.  The Limits of Resolution.  CRC Press, 2016.
  * [18] JE Fernandez Del Rio and MF Catedra-Perez.  The matrix pencil method for two-dimensional direction of arrival estimation employing an l-shaped array.  IEEE transactions on antennas and propagation, 45(11):1693–1694, 1997.
  * [19] Laurent Demanet and Nam Nguyen.  The recoverability limit for superresolution via sparsity.  arXiv preprint arXiv:1502.01385, 2015.
  * [20] Justin Demmerle, Eva Wegel, Lothar Schermelleh, and Ian M Dobbie.  Assessing resolution in super-resolution imaging.  Methods, 88:3–10, 2015.
  * [21] Arnold Jan Den Dekker and A Van den Bos.  Resolution: a survey.  JOSA A, 14(3):547–557, 1997.
  * [22] Quentin Denoyelle, Vincent Duval, and Gabriel Peyré.  Support recovery for sparse super-resolution of positive measures.  Journal of Fourier Analysis and Applications, 23(5):1153–1194, 2017\.
  * [23] G. Toraldo Di Francia.  Resolving power and information.  JOSA, 45(7):497–501, 1955.
  * [24] David L. Donoho.  Superresolution via sparsity constraints.  SIAM journal on mathematical analysis, 23(5):1309–1331, 1992.
  * [25] Vincent Duval and Gabriel Peyré.  Exact support recovery for sparse spikes deconvolution.  Foundations of Computational Mathematics, 15(5):1315–1355, 2015\.
  * [26] Walter Gautschi.  On inverses of vandermonde and confluent vandermonde matrices.  Numerische Mathematik, 4(1):117–123, 1962.
  * [27] Keyong Han and Arye Nehorai.  Improved source number detection and direction estimation with nested arrays and ulas using jackknifing.  IEEE Transactions on Signal Processing, 61(23):6118–6128, 2013\.
  * [28] Zhaoshui He, Andrzej Cichocki, Shengli Xie, and Kyuwan Choi.  Detecting the number of clusters in n-way probabilistic clustering.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 32(11):2006–2021, 2010.
  * [29] S.W. Hell and J. Wichmann.  Breaking the diffraction resolution limit by stimulated emission: stimulated-emission-depletion fluorescence microscopy.  Opt. Lett., 19:780–782, 1994.
  * [30] C Helstrom.  The detection and resolution of optical signals.  IEEE Transactions on Information Theory, 10(4):275–287, 1964.
  * [31] Carl W Helstrom.  Detection and resolution of incoherent objects by a background-limited optical system.  JOSA, 59(2):164–175, 1969.
  * [32] S.T. Hess, T.P.K. Girirajan, and M.D. Mason.  Ultra-high resolution imaging by fluorescence photoactivation localization microscopy.  Biophys. J., 91:4258–4272, 2006.
  * [33] William V Houston.  A compound interferometer for fine structure work.  Physical Review, 29(3):478, 1927.
  * [34] Richard L Johnson and Gina E Miner.  An operational system implementation of the ESPIRIT DF algorithm.  IEEE transactions on aerospace and electronic systems, 27(1):159–166, 1991.
  * [35] Shohei Kikuchi, Hiroyuki Tsuji, and Akira Sano.  Pair-matching method for estimating 2-d angle of arrival with a cross-correlation matrix.  IEEE Antennas and Wireless Propagation Letters, 5:35–40, 2006.
  * [36] DN Lawley.  Tests of significance for the latent roots of covariance and correlation matrices.  biometrika, 43(1/2):128–136, 1956.
  * [37] Weilin Li and Wenjing Liao.  Stable super-resolution limit and smallest singular value of restricted fourier matrices.  Applied and Computational Harmonic Analysis, 51:118–156, 2021.
  * [38] Weilin Li, Wenjing Liao, and Albert Fannjiang.  Super-resolution limit of the ESPRIT algorithm.  IEEE Transactions on Information Theory, 66(7):4593–4608, 2020\.
  * [39] Wenjing Liao.  MUSIC for multidimensional spectral estimation: stability and super-resolution.  IEEE transactions on signal processing, 63(23):6395–6406, 2015\.
  * [40] Wenjing Liao and Albert C. Fannjiang.  MUSIC for single-snapshot spectral estimation: Stability and super-resolution.  Applied and Computational Harmonic Analysis, 40(1):33–67, 2016\.
  * [41] Ping Liu.  Mathematical Theory of Computational Resolution Limit and Efficient Fast Algorithms for Super-Resolution.  Hong Kong University of Science and Technology (Hong Kong), 2021.
  * [42] Ping Liu and Hai Zhang.  A mathematical theory of computational resolution limit in multi-dimensional spaces.  Inverse Problems, 37(10):104001, 2021.
  * [43] Ping Liu and Hai Zhang.  A theory of computational resolution limit for line spectral estimation.  IEEE Transactions on Information Theory, 67(7):4812–4827, 2021\.
  * [44] Ping Liu and Hai Zhang.  A mathematical theory of computational resolution limit in one dimension.  Applied and Computational Harmonic Analysis, 56:402–446, 2022.
  * [45] Ping Liu and Hai Zhang.  A measurement decoupling based fast algorithm for super-resolving point sources with multi-cluster structure.  arXiv preprint arXiv:2204.00469, 2022.
  * [46] Tsung-Hsien Liu and Jerry M Mendel.  Azimuth and elevation direction finding using arbitrary array geometries.  IEEE Transactions on Signal Processing, 46(7):2061–2065, 1998.
  * [47] Leon B Lucy.  Resolution limits for deconvolved images.  The Astronomical Journal, 104:1260–1265, 1992.
  * [48] Leon B Lucy.  Statistical limits to super resolution.  Astronomy and Astrophysics, 261:706, 1992.
  * [49] Ankur Moitra.  Super-resolution, extremal functions and the condition number of vandermonde matrices.  In Proceedings of the Forty-seventh Annual ACM Symposium on Theory of Computing, STOC ’15, pages 821–830, New York, NY, USA, 2015. ACM.
  * [50] Veniamin I Morgenshtern.  Super-resolution of positive sources on an arbitrarily fine grid.  arXiv preprint arXiv:2005.06756, 2020.
  * [51] Veniamin I. Morgenshtern and Emmanuel J. Candès.  Super-resolution of positive sources: The discrete setup.  SIAM Journal on Imaging Sciences, 9(1):412–444, 2016.
  * [52] David W Pentico.  Assignment problems: A golden anniversary survey.  European Journal of Operational Research, 176(2):774–793, 2007\.
  * [53] Clarice Poon and Gabriel Peyré.  Multidimensional sparse super-resolution.  SIAM Journal on Mathematical Analysis, 51(1):1–44, 2019.
  * [54] Sripad Ram, E Sally Ward, and Raimund J Ober.  Beyond rayleigh’s criterion: a resolution measure with application to single-molecule microscopy.  Proceedings of the National Academy of Sciences, 103(12):4457–4462, 2006.
  * [55] Lord Rayleigh.  Xxxi. investigations in optics, with special reference to the spectroscope.  The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 8(49):261–274, 1879.
  * [56] Jorma Rissanen.  Modeling by shortest data description.  Automatica, 14(5):465–471, 1978.