# Spectral Properties of Infinitely Smooth Kernel Matrices in the Single
Cluster Limit, with Applications to Multivariate Super-Resolution

Nuha Diab  Department of Applied Mathematics, Tel Aviv University, Tel Aviv,
Israel  Dmitry Batenkov  Department of Applied Mathematics, Tel Aviv
University, Tel Aviv, Israel

###### Abstract

We study the spectral properties of infinitely smooth multivariate kernel
matrices when the nodes form a single cluster. We show that the geometry of
the nodes plays an important role in the scaling of the eigenvalues of these
kernel matrices. For the multivariate Dirichlet kernel matrix, we establish a
criterion for the sampling set ensuring precise scaling of eigenvalues.
Additionally, we identify specific sampling sets that satisfy this criterion.
Finally, we discuss the implications of these results for the problem of
super-resolution, i.e. stable recovery of sparse measures from bandlimited
Fourier measurements.

##  1 Introduction

Super-Resolution Problem. The problem of super-resolution (SR) is to recover
the fine details of a signal from inherently low resolution measurements in
the frequency domain. The SR problem has a variety of applications in signal
processing, imaging, optics, inverse scattering, and data analysis problems
[9, 27]. In the last years, a particular SR model has received considerable
attention, where the signal is modelled as a sparse measure [8]. Let
$\mathbb{T}^{d}:=[-\pi,\pi)^{d}\equiv\left(\mathbb{R}\mod 2\pi\right)^{d}$,
let $\delta(\cdot)$ denote the Dirac delta distribution, and consider the
signal

$f(t):=\sum_{j=1}^{n}\alpha_{j}\delta(t-x_{j}),\quad
t,x_{j}\in\mathbb{T}^{d},\quad\alpha_{j}\in\mathbb{C},$ (1)

which we aim to recover. We refer to $\\{x_{j}\\}$ as the ”nodes” and to
$\\{\alpha_{j}\\}$ as the ”coefficients” or “amplitudes”. The noisy
measurements are given by

$\hat{f}(\omega):=\sum_{j=1}^{n}\alpha_{j}e^{\imath\langle\omega,x_{j}\rangle}+\epsilon(\omega),\quad\omega\in{\cal
G}_{N}:=\\{-N+1,\dots,-1,0,1,\dots,N-1\\}^{d}\subset\mathbb{Z}^{d},\quad|\epsilon(\omega)|\leq\epsilon.$
(2)

Moreover, define $\text{SRF}:=\frac{1}{N\Delta}$, where $\Delta$ represents
the smallest distance (e.g. in the infinity norm) between the nodes. It is
well-known that, at least in the one-dimensional case, this quantity controls
the numerical stability of the problem. The super-resolution regime is when
$\text{SRF}\gg 1$; while the well-separated case is when $\text{SRF}\leq 1$.
In the super-resolution regime, nodes are organized in clusters, with a
distance of order $\Delta$ between any two nodes within a cluster.

One primary objective in theory of super-resolution is to define optimal
bounds for reconstruction errors, commonly referred to as min-max error
bounds, or the computational resolution limits [14, 28, 25, 12, 8]. These
bounds are achieved by the (non-tractable) most effective algorithm under the
worst-case scenario. Furthermore, these stability bounds serve to validate the
optimality of different tractable solution methods, offering assurances
regarding the methods’ performance [26, 20, 21]. It has been showed previously
that the min-max bounds are related to the smallest singular value of the
Vandermonde matrix $U$ of the system [12], which is defined as

$U(x_{1},\dots,x_{n};{\cal
G}_{N}):=\big{[}\exp{\imath\langle\omega,x_{j}\rangle}\big{]}_{\omega\in{\cal
G}_{N}}^{j=1,\dots,n}\in\mathbb{C}^{|{\cal G}_{N}|\times n}.$ (3)

Extensive research has been conducted on the one-dimensional scenario,
covering both well-separated and super-resolution regimes [7, 25, 11, 33]. Let
$\ell$ be the number of nodes in the largest cluster. Then, in the super-
resolution regime, the smallest singular value of the Vandermonde matrix
scales like $\text{SRF}^{1-\ell}$, resulting in “on-grid” min-max bounds to be
of order $\text{SRF}^{2\ell-1}\epsilon$ [5, 25], (cf. similar off-grid bounds
[8, 29]). In the multidimensional context, several results are available for
the well-separated regime as well as the super-resolution regime [28, 23, 22,
16, 34], under specific conditions related to the unknown nodes. A key
distinction between the one-dimensional and multidimensional cases lies in the
spectrum of the Vandermonde matrix, wherein the geometry of the nodes plays a
crucial role. As we shall show in this paper, and consistent with the
observations in [22], it is not solely the distance between the nodes that
determines the spectral properties of the Vandermonde matrix, at least in the
near-colliding limit, but also the algebraic variety on which these nodes are
situated.

Kernel matrices in the flat limit. Following the previous discussion, let us
consider the Gramian matrix $D_{N}:=\frac{1}{{(2N)}^{d}}U^{*}U$, which is the
“kernel matrix” for the multidimensional Dirichlet kernel. Kernel matrices are
of importance in various fields, including scattered data approximation and
machine learning, see e.g. [1, 41] and references therein. In the context of
super-resolution, if all the nodes $\mathcal{X}=\\{x_{j}\\}_{j=1}^{n}$ form a
single cluster and $\text{SRF}\gg 1$, then the kernel matrix $D_{N}$ can be
considered in the so-called “flat limit”, a term introduced in [15] in the
context of radial basis function (RBF) interpolation. The flat limit was
investigated recently in various publications, from among those [1] and [41]
being most relevant to our work. While [41] deals exclusively with RBF kernels
and derives asymptotic eigenvalue scaling, [1] extends those results to
arbitrary smooth kernels, albeit under the assumption of the points not lying
on any low-dimensional algebraic variety. The approach from [41] is primarily
based on Micchelli’s lemma [30], and we have adopted this approach previously
in the one-dimensional case [7].

Contributions. The main result of this paper, Theorem 1, derives the
asymptotic scaling of eigenvalues of a kernel matrix induced by an infinitely
smooth kernel in the flat limit, generalizing the corresponding result of [1]
to arbitrary geometry. As a corollary, we derive the asymptotic scaling of the
eigenvalues of $D_{N}$, leading to the appropriate scaling of the singular
values of the Vandermonde matrix (3), when the nodes are arranged within a
single cluster. In this context, a ”single cluster” is characterized by the
condition $N\Delta\ll 1$. We therefore obtain a generalization of the results
from [23] which considered mainly the “single line” geometry, as well as the
extension of [7, Theorem 2.3] to the multivariate case. See also Remark 6.

In more detail, suppose that $x_{j}=\Delta y_{j}$ with $\epsilon=N\Delta$, and
take $\epsilon\to 0$. We show that the number of eigenvalues of
$D_{N}(\\{x_{1},\dots,x_{n}\\})$ decaying like $\epsilon^{2k}$, is precisely
equal to a certain number $t_{k}$ which depends on the algebraic properties of
the original nodes $\\{y_{j}\\}$ (see Theorem 5). In particular, when all the
nodes belong to a one-dimensional affine subspace (a “single-line” geometry),
then $t_{k}=1$, which implies $\sigma_{\min}(U)\approx
N^{\frac{d}{2}}\epsilon^{n-1}$ (as derived in [23]). On the other hand, when
$\\{y_{j}\\}$ are in a general position, then, as long as $n\geq{k+d\choose
d}$, we have $t_{k}={k+d-1\choose d-1}$, where $d$ is the dimension of the
space. These numbers coincide for $d=1$, but may differ drastically for $d>1$.
For example, if $d=2$ then $\sigma_{\min}(U)$ may be as large as
$N\epsilon^{\mathcal{O}(\sqrt{n})}$. Our results are slightly more general and
provide the eigenvalue scaling for arbitrary (symmetric) sampling sets
$\omega\in\mathcal{S}$, with the upper bounds the same as described above;
however we can prove matching lower bounds only for certain $\mathcal{S}$
satisfying the so-called _Geometric Characterization Condition_ (in particular
satisfied by ${\cal G}_{N}$), see Definition 8 below.

To demonstrate the relevance of our findings to the multivariate SR problem,
in Section 6 we also report on a numerical study of reconstructing the signals
from noisy measurements using the nonlinear least squares (NLS) method and the
multidimensional ESPRIT method [36]. Our results indeed suggest that the
asymptotic conditioning of the problem may vary from the “worst-case” single-
line scenario (previously considered in [28]), to the “best-case” general
position scenario. We are therefore confident that our findings will enhance
the understanding of super-resolution stability in high-dimensional settings
and can contribute to the analysis of current multivariate SR recovery methods
such as [13, 32, 38, 39, 36]. Furthermore, we believe these results can serve
as an important step towards establishing min-max error bounds in the general
multidimensional scenario, when the nodes form multiple clusters (where some
kind of multidimensional “confluent” Vandermonde matrices may be required, cf.
[3, 2, 4, 6]).

Organization of the paper. In Section 2 we establish some notation and
definitions. In Section 3 we present our main result, Theorem 1, which
describes the asymptotic scaling of the eigenvalues of kernel matrices in the
flat limit. In Section 4 we specialize the results to the Dirichlet kernel
matrix, showing, in particular, that the bounds of Theorem 1 are tight. In
Section 5 we derive the corresponding scaling of the singular values of the
Vandermonde matrix, which requires additional technical arguments. Finally,
numerical results are presented in Section 6.

Acknowledgements. We thank Konstantin Usevich, Benedikt Diederichs, and Weilin
Li for helpful discussions. This work has been supported by the Israel Science
Foundation, grant 1793/20, by a collaborative grant from the Volkswagen
Foundation, and by the Ariane de Rothschild Women Doctoral Program.

##  2 Preliminaries

We adapt some definitions and notation from [1]. For a multi-index
$\boldsymbol{\alpha}=(\alpha_{1},\dots,\alpha_{d})\in\mathbb{Z}^{d}_{+}$,
denote

$\big{|}\boldsymbol{\alpha}\big{|}:=\sum_{j=1}^{d}\alpha_{j},\quad\boldsymbol{\alpha}!:=\alpha_{1}!\cdot\cdot\cdot\alpha_{d}!,$

and let for each $j=0,1,2,\dots$

$\mathbb{P}_{j}:=\\{\boldsymbol{\alpha}\in\mathbb{Z}^{d}_{+}:\big{|}\boldsymbol{\alpha}\big{|}\leq
j\\},\quad\mathbb{H}_{j}:=\\{\boldsymbol{\alpha}\in\mathbb{Z}^{d}_{+}:\big{|}\boldsymbol{\alpha}\big{|}=j\\}=\mathbb{P}_{j}\setminus\mathbb{P}_{j-1},\quad\mathbb{P}_{-1}:=\emptyset.$

The cardinalities of these sets are given by:

$p_{j}:=\\#\mathbb{P}_{j}={j+d\choose d},\quad
h_{j}:=\\#\mathbb{H}_{j}={j+d-1\choose
d-1}=\\#\mathbb{P}_{j}-\\#\mathbb{P}_{j-1}.$ (4)

For $\boldsymbol{x}\in\mathbb{R}^{d}$ and
$\boldsymbol{\alpha}\in\mathbb{Z}^{d}_{+}$, the monomial
$\boldsymbol{x}^{\boldsymbol{\alpha}}$ is defined as
$\boldsymbol{x}^{\boldsymbol{\alpha}}:=x_{1}^{\alpha_{1}}\cdot\cdot\cdot
x_{d}^{\alpha_{d}}$, and its (total) degree is $|\boldsymbol{\alpha}|$. Let
$\prec$ be a fixed graded order on the set of multi-indices
$\mathbb{Z}^{d}_{+}$ (or, equivalently, on the set of monomials), for example
a graded lexicographic order or a graded reverse lexicographic order. It will
always be implied in the following.

######  Definition 1.

Let $\Omega\subseteq\mathbb{R}^{d}$ be a fixed open set. A kernel function
$\mathcal{K}:\Omega\times\Omega\rightarrow\mathbb{R}$ is said to belong to the
class $\mathcal{C}^{(r,r)}(\Omega)$ (or
$\mathcal{C}^{\mathbb{P}_{r}\times\mathbb{P}_{r}}(\Omega)$) when the partial
derivatives

$\mathcal{K}^{(\boldsymbol{\alpha,\beta})}(\boldsymbol{x},\boldsymbol{y}):=\frac{\partial^{|\boldsymbol{\alpha}|+|\boldsymbol{\beta}|}}{\partial
x_{1}^{\alpha_{1}}\cdot\cdot\cdot\partial x_{d}^{\alpha_{d}}\partial
y_{1}^{\beta_{1}}\cdot\cdot\cdot\partial
y_{d}^{\beta_{d}}}\mathcal{K}(\boldsymbol{x},\boldsymbol{y})$

exist and are continuous on $\Omega\times\Omega$, for all
$\boldsymbol{\alpha,\beta}\in\mathbb{P}_{r}$. When
$\mathcal{K}\in\mathcal{C}^{(r,r)}(\Omega)$ for all $r\in\mathbb{N}$, we write
$\mathcal{K}\in\mathcal{C}^{(\infty,\infty)}(\Omega)$.

######  Definition 2.

For a kernel $\mathcal{K}\in\mathcal{C}^{(r,r)}(\Omega)$, a finite set ${\cal
X}\subset\Omega$ and $\epsilon>0$ we define the scaled kernel
$\mathcal{K}_{\epsilon}(\boldsymbol{x},\boldsymbol{y}):=\mathcal{K}(\epsilon\boldsymbol{x},\epsilon\boldsymbol{y})$
and the kernel matrix

$K_{\epsilon}({\cal
X}):=\big{[}\mathcal{K}_{\epsilon}(\boldsymbol{x},\boldsymbol{x^{\prime}})\big{]}_{\boldsymbol{x},\boldsymbol{x^{\prime}}\in{\cal
X}}\in\mathbb{R}^{|{\cal X}|\times|{\cal X}|}.$ (5)

######  Definition 3 (Multivariate Vandermonde matrix).

For a finite set ${\cal
X}=\\{\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{n}\\}\subset\mathbb{R}^{d}$ of
nodes, and an $\prec$-ordered set of multi-indices
$\mathcal{A}=\\{\boldsymbol{\alpha}_{1},\dots,\boldsymbol{\alpha}_{m}\\}\subset\mathbb{Z}^{d}_{+}$
(i.e. $\boldsymbol{\alpha}_{1}\prec\dots\prec\boldsymbol{\alpha}_{m}$), we
define the multivariate Vandermonde matrix as

$V_{\mathcal{A}}=V_{\mathcal{A}}({\cal
X})=\big{[}(\boldsymbol{x}_{i})^{\boldsymbol{\alpha}_{j}}\big{]}_{1\leq i\leq
n}^{{1\leq j\leq m}}\in\mathbb{R}^{n\times m}.$ (6)

Notice that $V$ is a real Vandermonde matrix, as opposed to the (transpose of
the) complex Vandermonde matrix $U$ defined in (3). Let us further denote
$V_{\leq k}:=V_{\mathbb{P}_{k}}({\cal X})$ and
$V_{k}:=V_{\mathbb{H}_{k}}({\cal X})$. The matrix $V_{\leq
k}\in\mathbb{R}^{n\times p_{k}}$ can be partitioned into $k+1$ matrices
$V_{k}\in\mathbb{R}^{n\times h_{k}}$ arranged by increasing degree:

$V_{\leq k}=\big{[}V_{0}\quad V_{1}\quad\dots\quad V_{k}\big{]}.$ (7)

######  Definition 4.

Let $\mathcal{K}\in\mathcal{C}^{(r,r)}(\Omega)$, with
$\Omega\subseteq\mathbb{R}^{d}$ and $\mathbf{0}\in\Omega$. Let
$\mathcal{A},\mathcal{B}\subset\mathbb{Z}^{d}_{+}$ be two sets of multi-
indices satisfying $|\alpha|,|\beta|\leq r$ for all
$\alpha\in\mathcal{A},\beta\in\mathcal{B}$. The Wronskian matrix of
$\mathcal{K}$ is defined as:

$W^{\mathcal{K}}_{\mathcal{A},\mathcal{B}}=\bigg{[}{{\mathcal{K}^{(\alpha,\beta)}(0,0)}\over{\alpha!\beta!}}\bigg{]}_{\alpha\in\mathcal{A},\beta\in\mathcal{B}},$
(8)

where the rows and columns are indexed by multi-indices in $\mathcal{A}$ and
$\mathcal{B}$ according to the chosen ordering. In addition, denote by
$W^{\mathcal{K}}_{\leq
k}=W^{\mathcal{K}}_{\mathbb{P}_{k}}=W^{\mathcal{K}}_{\mathbb{P}_{k},\mathbb{P}_{k}}$.

The final piece of notation is central to our discussion, describing a
geometric property of ${\cal X}$.

######  Definition 5 (Discrete moment order, [40]).

The _discrete moment order_ of a set ${\cal X}$ consisting of $n$ pairwise
distinct points, is the smallest number $m=\mu(\mathcal{X})$ such that

$\ker(V_{\leq m}({\cal X})^{T})=\emptyset\iff\operatorname{rank}(V_{\leq
m}({\cal X}))=n.$

Note that such $m$ exists (is finite). Indeed, if ${\cal
X}=\\{x_{1},\dots,x_{n}\\}$, then $\mu({\cal X})\leq n-1$ because the
following Lagrange interpolation polynomials, each of total degree $n-1$, are
linearly independent:

$\biggl{\\{}\ell_{i}(x)=\prod_{j\neq
i}\frac{H_{j}(x)}{H_{j}(x_{i})}\biggr{\\}}_{i=1}^{n},$

where $H_{j}$ is some affine hyperplane containing $x_{j}$ and not containing
every other $x_{k},k\neq j$. Cf. [40, Definition 3],[17] and [41].

##  3 Eigenvalues of smooth kernel matrices in the flat limit

In this section, we consider the kernel matrix
$K_{\epsilon}=K_{\epsilon}({\cal X})$ for a fixed infinitely smooth kernel
$\mathcal{K}$ and a fixed set ${\cal X}\subset\mathbb{R}^{d}$ of $n$ pairwise
distinct points. Our main interest is in the asymptotic order of the decay of
the eigenvalues of $K_{\epsilon}$ in the flat limit $\epsilon\to 0$. We employ
the standard $O(\cdot)$ notation.

Our main result gives an upper bound on the decay rates of all the
eigenvalues.

######  Theorem 1.

Let $m=\mu({\cal X})$ be the discrete moment order of the set ${\cal X}$. For
any infinitely smooth kernel
$\mathcal{K}\in\mathcal{C^{(\infty,\infty)}}(\Omega)$ with
$\mathbf{0}\in\Omega$, and small enough $\epsilon$, the eigenvalues of
$K_{\epsilon}({\cal X})$ split into $m+1$ groups

$\displaystyle\lambda_{0,0}$
$\displaystyle=O(1),\quad\\{\lambda_{1,j}\\}_{j=1}^{t_{1}}=O(\epsilon^{2}),\quad\dots\quad,\quad\\{\lambda_{m,j}\\}_{j=1}^{t_{m}}=O(\epsilon^{2m}),$
$\displaystyle t_{k}$ $\displaystyle:=\operatorname{rank}(V_{\leq k}({\cal
X}))-\operatorname{rank}(V_{\leq k-1}({\cal X})),\;k=1,2,\dots,m.$

######  Remark 1.

Theorem (1) is valid also for kernels in $\mathcal{C}^{(m+1,m+1)}$.

###### Proof.

We start by writing the Taylor expansion of $K_{\epsilon}$ around 0 as
equation (55) in [1]:

$K_{\epsilon}=V_{\leq m}\Delta_{m}W\Delta_{m}V_{\leq
m}^{T}+\epsilon^{m+1}(V_{\leq
m}\Delta_{m}W_{1}(\epsilon)+W_{2}(\epsilon)\Delta_{m}V_{\leq
m}^{T})+\epsilon^{2{(m+1)}}W_{3}(\epsilon)$ (9)

where $W_{i}(\epsilon)=O(1)$, $W=W^{\mathcal{K}}_{\leq m}$ and
$\Delta_{m}(\epsilon)=\operatorname{diag}(1,\epsilon
I_{h_{1}},\dots,\epsilon^{m}I_{h_{m}})\in\mathbb{R}^{p_{m}\times p_{m}}$. For
the full derivation of (9), see section A in the Appendix.

Let $V=V_{\leq m}$ where $V=[V_{0}\dots V_{m}]=[V_{\leq m-1}\quad V_{m}]$, and
$t_{k}:=\operatorname{rank}(V_{\leq k})-\operatorname{rank}(V_{\leq k-1})$ as
defined in the theorem. In [1] only the generic case $t_{k}=h_{k}-h_{k-1}$ is
considered. Here we generalize Theorem 6.1 in [1] to the case where the
multivariate Vandermonde matrices $V_{\leq k}$ are not necessarily full rank.

Consider the full QR decomposition of $V_{\leq m-1}$ with column pivoting, as
elaborated in Section 5.4.2, specifically, formula (5.4.6) in the book [18]:

$V_{\leq m-1}P_{m-1}=[Q_{\leq m-1}\quad
Q_{\perp}]\begin{bmatrix}R_{11}&R_{12}\\\
0&0\end{bmatrix}=Q\begin{bmatrix}R\\\ 0\end{bmatrix}$

where $Q_{\leq m-1}\in\mathbb{R}^{n\times r_{m-1}}$,
$Q_{\perp}\in\mathbb{R}^{n\times t_{m}}$, $R_{11}\in\mathbb{R}^{r_{m-1}\times
r_{m-1}}$, $R\in\mathbb{R}^{r_{m-1}\times p_{m-1}}$, $r_{m-1}+t_{m}=n$ and
$P_{m-1}$ is a permutation matrix.
We have $[Q_{\leq m-1}\quad Q_{\perp}]^{T}V_{\leq m}P=\begin{bmatrix}R&Q_{\leq
m-1}^{T}V_{m}\\\ 0&Q_{\perp}^{T}V_{m}\end{bmatrix}$ and $Q_{\perp}^{T}V_{\leq
m}P=[0\quad Q_{\perp}^{T}V_{m}]$ where $P:=\begin{bmatrix}P_{m-1}&0\\\
0&I_{h_{m}}\end{bmatrix}$. Note that for every permutation matrix, we have
$P^{-1}=P^{T}$.

Thus we can rewrite the expansion (9) as follows:

$\displaystyle K_{\epsilon}$ $\displaystyle=\quad V_{\leq
m}PP^{T}\Delta_{m}PP^{T}WPP^{T}\Delta_{m}PP^{T}V_{\leq m}^{T}$
$\displaystyle\quad+\epsilon^{m+1}(V_{\leq
m}PP^{T}\Delta_{m}PP^{T}W_{1}(\epsilon)+W_{2}(\epsilon)PP^{T}\Delta_{m}PP^{T}V_{\leq
m}^{T})$ $\displaystyle\quad+\epsilon^{2{(m+1)}}W_{3}(\epsilon)$
$\displaystyle K_{\epsilon}$ $\displaystyle=\tilde{V}_{\leq
m}\tilde{\Delta}_{m}\tilde{W}\tilde{\Delta}_{m}\tilde{V}_{\leq
m}^{T}+\epsilon^{m+1}(\tilde{V}_{\leq
m}\tilde{\Delta}_{m}\tilde{W}_{1}(\epsilon)+\tilde{W}_{2}(\epsilon)\tilde{\Delta}_{m}\tilde{V}_{\leq
m}^{T})+\epsilon^{2{(m+1)}}W_{3}(\epsilon)$

where $\tilde{V}_{\leq m}:=V_{\leq m}P$ ,
$\tilde{\Delta}_{m}:=P^{T}\Delta_{m}P=\operatorname{diag}(\epsilon^{i_{1}}I_{h_{i_{1}}},...,\epsilon^{i_{m}}I_{h_{i_{m}}},\epsilon^{m}I_{h_{m}})$
due to Lemma 7, $\tilde{W}:=P^{T}WP$,
$\tilde{W}_{1}(\epsilon):=P^{T}W_{1}(\epsilon)$ and
$\tilde{W}_{2}(\epsilon):=W_{2}(\epsilon)P$.
Multiplying by $Q_{\perp}^{T}$ and its transpose we get:

$\displaystyle Q_{\perp}^{T}K_{\epsilon}Q_{\perp}$ $\displaystyle=\quad[0\quad
Q_{\perp}^{T}V_{m}]\tilde{\Delta}\tilde{W}\tilde{\Delta}[0\quad
Q_{\perp}^{T}V_{m}]^{T}$ $\displaystyle\quad+\epsilon^{m+1}[0\quad
Q_{\perp}^{T}V_{m}]\tilde{\Delta}\tilde{W}_{1}Q_{\perp}$
$\displaystyle\quad+\epsilon^{m+1}Q_{\perp}^{T}\tilde{W}_{2}\tilde{\Delta}[0\quad
Q_{\perp}^{T}V_{m}]^{T}$
$\displaystyle\quad+(\epsilon^{m+1}Q_{\perp}^{T}W_{3}Q_{\perp}\epsilon^{m+1}).$

We can also write the above as:

$\displaystyle Q_{\perp}^{T}K_{\epsilon}Q_{\perp}$ $\displaystyle=\quad
O(\epsilon^{2m})Q_{\perp}^{T}V_{m}\hat{W}V_{m}^{T}Q_{\perp}$
$\displaystyle\quad+O(\epsilon^{2m+1})(Q_{\perp}^{T}V_{m}\hat{W_{1}}Q_{\perp}$
$\displaystyle\quad+Q^{T}_{\perp}\hat{W_{2}}V_{m}Q_{\perp})$
$\displaystyle\quad+O(\epsilon^{2{(m+1)}}),$

where $\hat{W}$ is the lower right $h_{m}\times h_{m}$ sub-matrix of
$\tilde{W}$. Thus, by the Courant-Fischer principle, $K_{\epsilon}$ has at
least $t_{m}=\dim(Q_{\perp})$ eigenvalues which decay as fast as
$O(\epsilon^{2m})$.

For any $k=1,...,m-1$, let $V_{\leq k-1}P_{k-1}=[Q_{\leq k-1}\quad
Q_{\perp,k-1}]\begin{bmatrix}R_{\leq k-1}\\\ 0\end{bmatrix}$ be the full QR
decomposition of $V_{\leq k}$ with pivoting, where $Q_{\leq
k-1}\in\mathbb{R}^{n\times[n-(t_{k}+\dots+t_{m})]}$,
$Q_{\perp,k-1}\in\mathbb{R}^{n\times[t_{k}+\dots+t_{m}]}$ and $R_{\leq
k-1}\in\mathbb{R}^{[n-(t_{k}+\dots+t_{m})]\times p_{k-1}}$.

We can apply the same argument on $Q_{\perp,k-1}$ and $V_{\leq k-1}$ by
induction, then we get that there are at least $t_{k}+...+t_{m}$ eigenvalues
which decay as fast as $O(\epsilon^{2k})$ and there are at least
$t_{k+1}+...+t_{m}$ eigenvalues which decay as fast as
$O(\epsilon^{2({k+1})})$. Therefore, there are at least $t_{k}$ eigenvalues
which decay as fast as $O(\epsilon^{2{k}})$. Since $\sum_{i=0}^{m}t_{i}=n$
(assume $t_{0}=0$) we get that there are exactly $t_{k}$ eigenvalues decaying
(at least) like $\epsilon^{2k}$ when $\epsilon$ goes to 0. ∎

The bounds of Theorem 1 are not necessarily tight. As we show in the next
lemma, the exactness of the asymptotic orders depends on the analytic
properties of the kernel $\mathcal{K}$ and the geometry of ${\cal X}$,
therefore additional assumptions ought to be made on both of those objects in
order to obtain matching lower bounds.

######  Lemma 1.

For $\mathcal{K}\in\mathcal{C^{(\infty,\infty)}}(\Omega)$ and ${\cal
X}\subset\mathbb{R}^{d}$, there exists $C=C({\cal
X},\mathcal{K})\in\mathbb{R}$ such that

$\det(K_{\epsilon}({\cal
X}))=\epsilon^{2\sum_{i=0}^{m}it_{i}}(C+O(\epsilon)),\quad\epsilon\to 0.$

In particular, $C=\det\left(\tilde{R}W\tilde{R}^{T}\right)$ where
$\tilde{R}\in\mathbb{R}^{n\times p_{m}}$ is full rank, and
$W=W^{\mathcal{K}}_{\leq m}\in\mathbb{R}^{p_{m}\times p_{m}}$ is the Wronskian
matrix of $\mathcal{K}$ as in Definition 4.

###### Proof.

Consider the matrices $Q_{\perp,k-1}$ defined in the proof of Theorem 1. Let
$\mathcal{M}_{k}\subseteq range(Q_{\perp,k-1})$ and $\mathcal{M}_{k}\perp
range(Q_{\perp,k})$ where $\mathcal{M}_{k}\in\mathbb{R}^{n\times t_{k}}$ and
$Q_{\perp,-1}=\mathbb{R}^{n\times n}$. Define

$\mathcal{Q}_{\perp}:=\big{(}\mathcal{M}_{0},\mathcal{M}_{1},...,\mathcal{M}_{m}\big{)}\in\mathbb{R}^{n\times
n}$

It is easy to see that $\mathcal{Q}_{\perp}^{T}V_{\leq m}$ is upper
triangluar, with the block-diagonal part given by

$\tilde{R}:=\operatorname{blkdiag}\left(\mathcal{Q}_{\perp}^{T}V_{\leq
m}\right)=\operatorname{diag}\bigl{\\{}\mathcal{M}_{0}^{T}V_{0},...,\mathcal{M}_{m}^{T}V_{m}\bigr{\\}}.$

Each diagonal block of $\tilde{R}$ is full rank by definition, therefore
$\tilde{R}$ is full rank.

Therefore by [1, Lemma 6.4] we have

$E_{n}^{-1}\mathcal{Q}_{\perp}^{T}V_{\leq m}\Delta_{m}=\tilde{R}+O(\epsilon)$

where $E_{n}=\operatorname{diag}\\{1,\epsilon
I_{t_{1}},...,\epsilon^{m}I_{t_{m}}\\}\in\mathbb{R}^{n\times n}$ and
$\tilde{R}\in\mathbb{R}^{n\times p_{m}}$.
Thus, using (9), we get that

$\displaystyle
E_{n}^{-1}\mathcal{Q}_{\perp}^{T}K_{\epsilon}\mathcal{Q}_{\perp}E_{n}^{-1}$
$\displaystyle=(\tilde{R}+O(\epsilon))W(\tilde{R}^{T}+O(\epsilon))$
$\displaystyle+\epsilon^{m+1}((\tilde{R}+O(\epsilon))W_{1}(\epsilon)\mathcal{Q}_{\perp}E_{n}^{-1}$
$\displaystyle+E_{n}^{-1}\mathcal{Q}_{\perp}^{T}W_{2}(\epsilon)(\tilde{R}^{T}+O(\epsilon)))$
$\displaystyle+\epsilon^{2{(m+1)}}E_{n}^{-1}\mathcal{Q}_{\perp}^{T}W_{3}(\epsilon)\mathcal{Q}_{\perp}E_{n}^{-1}$
$\displaystyle=\tilde{R}W\tilde{R}^{T}+O(\epsilon)$

where the last equality follows from
$\epsilon^{2{(m+1)}}E_{n}^{-1}=O(\epsilon)$ and
$\tilde{R}W\tilde{R}^{T}\in\mathbb{R}^{n\times n}$ . This implies that
$\epsilon^{-2\sum_{j=0}^{m}jt_{j}}\det(K_{\epsilon})=\det(\tilde{R}W\tilde{R}^{T}+O(\epsilon))=\det(\tilde{R}W\tilde{R}^{T})+O(\epsilon)$
(by Lemma (8) in the Appendix) with $C=\det(\tilde{R}W\tilde{R}^{T})$.

∎

######  Corollary 1.

Suppose $\mathcal{K}$ is symmetric and analytic (in both variables) in the
neighborhood of $(\mathbf{0},\mathbf{0})$. The scaling of the eigenvalues as
given in Theorem 1 is exact, if and only if $\det(\tilde{R}W\tilde{R}^{T})\neq
0$.

###### Proof.

By Theorem 1 and [1, Theorem 2.9], the eigenvalues of $K_{\epsilon}({\cal X})$
satisfy

$\lambda_{0,0}=\epsilon^{0}(\tilde{\lambda}_{0,0}+O(\epsilon)),\quad\\{\lambda_{1,j}\\}_{j=1}^{t_{1}}=\\{\epsilon^{2}(\tilde{\lambda}_{1,j}+O(\epsilon))\\}_{j=1}^{t_{1}},...,\\{\lambda_{m,j}\\}_{j=1}^{t_{m}}=\\{\epsilon^{2m}(\tilde{\lambda}_{m,j}+O(\epsilon))\\}_{j=1}^{t_{m}},$

where $t_{k}:=\operatorname{rank}(V_{\leq k})-\operatorname{rank}(V_{\leq
k-1})$ and $m=\mu({\cal X})$, and $\tilde{\lambda}_{s,j}$ do not depend on
$\epsilon$ (but may be zero). By Lemma 1 we have that

$\big{|}K_{\epsilon}({\cal
X})\big{|}=\epsilon^{2\sum_{i=0}^{m}it_{i}}\biggl{\\{}\big{|}\tilde{R}W\tilde{R}^{T}\big{|}+O(\epsilon)\biggr{\\}}.$

On the other hand, the determinant is the product of eigenvalues, therefore

$\displaystyle\big{|}K_{\epsilon}({\cal X})\big{|}$
$\displaystyle=\prod_{j=1}^{n}\lambda_{j}=\lambda_{0,0}\prod_{j=1}^{t_{1}}\lambda_{1,j}...\prod_{j=1}^{t_{m}}\lambda_{m,j}=(\tilde{\lambda}_{0,0}+O(\epsilon))\prod_{j=1}^{t_{1}}\epsilon^{2}(\tilde{\lambda}_{1,j}+O(\epsilon))...\prod_{j=1}^{t_{m}}\epsilon^{2m}(\tilde{\lambda}_{m,j}+O(\epsilon))$
$\displaystyle=\epsilon^{2\sum_{i=0}^{m}it_{i}}\biggl{\\{}\prod_{s=1}^{m}\prod_{j=1}^{t_{s}}\tilde{\lambda}_{s,j}+O(\epsilon)\biggr{\\}}.$

Now we just compare the two expressions above. In one direction, if the
scaling is exact then $\tilde{\lambda}_{s,j}\neq 0$ for all $0\leq s\leq m$
and $0\leq j\leq t_{s}$, which implies
$\left|\tilde{R}W\tilde{R}^{T}\right|\neq 0$. In the other direction, if there
are $\ell\geq 1$ zeros
$\tilde{\lambda}_{k}\in\\{\tilde{\lambda}_{s,j}\\}_{s=1,\dots,m}^{j=1,\dots,t_{m}}$
and $\tilde{\lambda}_{k}=0$ for $k=1,...,\ell$, then
$\big{|}K_{\epsilon}\big{|}=\epsilon^{\ell+2\sum_{i=0}^{m}it_{i}}(\prod_{k=\ell+1}^{n}\tilde{\lambda}_{k}+O(\epsilon))$,
forcing $\left|\tilde{R}W\tilde{R}^{T}\right|=0$. ∎

##  4 Eigenvalues of the multidimensional Dirichlet kernel

In this section we consider the specific case of the multidimensional
Dirichlet kernel, which is directly motivated by the problem of super-
resolution as outlined in the Introduction. Let ${\cal
S}\subset\mathbb{R}^{d}$ denote a finite set of frequencies, symmetric around
the origin, where the spectral data (2) is available. We first derive a
generic sufficient condition on ${\cal S}$ which ensures that the scaling
order of Theorem 1 is exact. Later we show that this condition is satisfied by
${\cal S}={\cal G}_{N}$, the regular symmetric grid in $\mathbb{R}^{d}$ of
side length $2N+1$, for large enough $N$.

######  Definition 6.

For ${\cal S}\subset\mathbb{R}^{d}$, symmetric around
$\mathbf{0}\in\mathbb{R}^{d}$, the (real-valued) Dirichlet kernel
$\mathcal{D}_{{\cal S}}$ is defined as follows:

$\mathcal{D}_{{\cal S}}(x)=\sum_{\omega\in{\cal S}}\exp\left(\imath\langle
x,\omega\rangle\right).$

By a slight abuse of notation, we also denote by $\mathcal{D}_{{\cal S}}$ the
corresponding bi-multivariate kernel:

$\mathcal{D}_{{\cal S}}(x,y)=\mathcal{D}_{{\cal S}}(x-y)=\sum_{\omega\in{\cal
S}}\exp\left(\imath\langle x-y,\omega\rangle\right).$ (10)

Given ${\cal X}\subset\mathbb{R}^{d}$ and $\epsilon>0$, the corresponding
Dirichlet kernel matrix and its flat limit version are denoted as:

$\displaystyle\begin{split}D({\cal X},{\cal S})&:=\big{[}\mathcal{D}_{{\cal
S}}(x,x^{\prime})\big{]}_{x,x^{\prime}\in{\cal X}},\\\ D_{\epsilon}({\cal
X},{\cal S})&:=D(\epsilon{\cal X},{\cal S})=\big{[}\mathcal{D}_{{\cal
S}}(\epsilon x,\epsilon x^{\prime})\big{]}_{x,x^{\prime}\in{\cal
X}}.\end{split}$ (11)

In addition, if ${\cal S}={\cal S}_{1}\times...\times{\cal S}_{d}$ is a tensor
product grid, we have $\mathcal{D}_{{\cal
S}}(x,y)=\prod_{j=1}^{d}\mathcal{D}_{{\cal S}_{j}}(x_{j},y_{j})$.

###  4.1 A nondegeneracy condition

Note that $\mathcal{D}_{{\cal S}}$ is a translation invariant, symmetric and
analytic kernel. Recall the definition (8) of the Wronskian matrix for an
arbitrary kernel. Set $W:=W^{\mathcal{D_{S}}}_{\leq m}$. By Corollary 1, in
order to show tightness, it is sufficient to demonstrate that
$C:=\det(\tilde{R}W\tilde{R}^{T})\neq 0$. In the following lemma, we show that
this holds under a certain (more or less natural) non-degeneracy assumption on
the set ${\cal S}$.

######  Lemma 2.

For $W:=W_{\leq m}$ as defined above, $\det(\tilde{R}W\tilde{R}^{T})>0$ holds
whenever $\operatorname{rank}(V_{\leq m}({\cal S}))\geq n$.

######  Remark 2.

Note that contrary to Section 3, the Vandermonde matrix $V_{\leq m}$ is here
evaluated on the set ${\cal S}$ (and not on ${\cal X}$), while the polynomial
total degree is still $m=\mu({\cal X})$.

###### Proof of Lemma 2.

We first show that $W$ is positive semidefinite. Directly from (8) and (10),
we have

$\displaystyle W$ $\displaystyle=\bigg{[}{{\mathcal{D}_{{\cal
S}}^{(\alpha,\beta)}(0,0)}\over{\alpha!\beta!}}\bigg{]}_{\alpha,\beta\in\mathbb{P}_{m}}=\bigg{[}\sum_{\omega\in{\cal
S}}{{(i\omega)^{\alpha}(-i\omega)^{\beta}}\over{\alpha!\beta!}}\bigg{]}_{\alpha,\beta\in\mathbb{P}_{m}}$
$\displaystyle=\bigg{[}\sum_{\omega=(\omega_{1},\dots,\omega_{d})\in{\cal
S}}{{(i)^{|\alpha|}(-i)^{|\beta|}\omega_{1}^{\alpha_{1}+\beta_{1}}\cdot...\cdot\omega_{d}^{\alpha_{d}+\beta_{d}}}\over{\alpha!\beta!}}\bigg{]}_{\alpha,\beta\in\mathbb{P}_{m}}$
$\displaystyle=FV^{T}_{\leq m}({\cal S})V_{\leq m}({\cal S})F^{*},$

where $V_{\leq m}({\cal S})=\big{[}{(\omega)}^{\alpha}\big{]}_{\omega\in{\cal
S}}^{\alpha\in\mathbb{P}_{m}}$ is the multivariate Vandermonde matrix and
$F:=\operatorname{diag}\left(\frac{i^{|\alpha|}}{\alpha!}\right)_{\alpha\in\mathbb{P}_{m}}$.
Let $B:=FV^{T}_{\leq m}({\cal S})$, and consider $B=U\Sigma T^{*}$ the
singular value decomposition of $B$. We therefore have
$W=BB^{*}=U\Sigma^{2}U^{*}$, implying that $W$ is positive semidefinite (since
$W$ is Hermitian and has non-negative eigenvalues). Since $\tilde{R}$ is full
rank by Lemma 1, it follows that $\tilde{R}W\tilde{R}^{T}$ is positive
semidefinite as well, implying that $\det(\tilde{R}W\tilde{R}^{T})\geq 0$. To
ensure $\det(\tilde{R}W\tilde{R}^{T})>0$, we require that
$\tilde{R}W\tilde{R}^{T}$ be positive definite. The latter condition is
satisfied when $\operatorname{rank}(\tilde{R}W\tilde{R}^{T})\geq n$, which is
equivalent to $\operatorname{rank}(W)=\operatorname{rank}(V_{\leq m}({\cal
S}))\geq n$. ∎

######  Remark 3.

For $d=1$, we have that $m=n-1$ and $p_{n-1}={n\choose 1}=n$, meaning that the
condition $\operatorname{rank}(V_{\leq n-1}({\cal S}))\geq n$ in Lemma (2) is
satisfied as long as the sampling set ${\cal S}$ contains at least $n$
distinct points.

######  Remark 4.

Given the condition $\operatorname{rank}(V_{\leq m}({\cal S}))\geq n$ in Lemma
(2), we aim to find a sampling set ${\cal S}_{*}$ with $p_{\ell}$ points where
$p_{\ell-1}<n$ and $p_{\ell}\geq n$ such that
$\operatorname{rank}(V_{\leq\ell}({\cal S}_{*}))\geq n$ which is independent
of $m=\mu({\cal X})$ (the geometry of the nodes).

By applying Corollary 1, we immediately obtain the following result.

######  Lemma 3.

Assume that the sampling set ${\cal S}$ satisfies the condition
$\operatorname{rank}(V_{\leq m}({\cal S}))\geq n$. Then the eigenvalues of
$D_{\epsilon}({\cal X},{\cal S})$ split into $m+1$ groups

$\lambda_{0,0}=\epsilon^{0}(\tilde{\lambda}_{0,j}+O(\epsilon)),\quad\\{\lambda_{1,j}\\}_{j=1}^{t_{1}}=\\{\epsilon^{2}(\tilde{\lambda}_{1,j}+O(\epsilon))\\}_{j=1}^{t_{1}},...,\\{\lambda_{m,j}\\}_{j=1}^{t_{m}}=\\{\epsilon^{2m}(\tilde{\lambda}_{m,j}+O(\epsilon))\\}_{j=1}^{t_{m}},$

where $t_{k}:=\operatorname{rank}(V_{\leq k}({\cal
X}))-\operatorname{rank}(V_{\leq k-1}({\cal X}))$ and
$\tilde{\lambda}_{s,j}\neq 0$ for all $0\leq s\leq m$ and $0\leq j\leq t_{s}$.

###  4.2 Geometric characterization condition and the uniform grid

The rank condition in Lemma 2 is nontrivial to verify for a given set ${\cal
S}$. In this section we provide an explicit example of such a set, namely, the
uniform grid ${\cal G}_{N}$ in $\mathbb{R}^{d}$. Our approach is a
straightforward application of well-known results in multivariate polynomial
interpolation, adapted to our setting.

We first recall some definitions and results from [37].

######  Definition 7 ([37], p.194).

Let $\mathcal{P}$ be a linear subspace of the polynomial ring
$\Pi:=\mathbb{R}\left[x_{1},...,x_{d}\right]$. The polynomial interpolation
problem with respect to a (finite) set of linearly independent functionals
$\Theta\subset\Pi^{\prime}$ ($\Pi^{\prime}$ is the dual space) is said to be
poised for $\mathcal{P}$, if for any $Y\in\mathbb{R}^{\\#\Theta}$ there exists
a unique $f\in\mathcal{P}$ such that $\Theta f=Y$.

######  Theorem 2 ([37] p.194).

For $\mathcal{P}\subset\Pi$, and a finite set $\Theta\subset\Pi^{\prime}$ as
above, the following statements are equivalent:

  1. 1.

The polynomial interpolation problem with respect to $\Theta$ is poised for
$\mathcal{P}$.

  2. 2.

$\dim\mathcal{P}=\\#\Theta$ and $\Theta P=\big{[}\theta p:\theta\in\Theta,p\in
P\big{]}$ satisfies $|\Theta P|\neq 0$ for any basis $P$ of $\mathcal{P}$.

For $\Theta_{{\cal S}}:=\\{\delta_{\xi}\\}_{\xi\in{\cal S}}$ the set of point
evaluation functionals on ${\cal S}$, and the standard monomial basis $P_{n}$
of polynomials of total degree $\leq n$, we have in particular $\dim
P_{n}=\\#{\cal S}={{n+d}\choose d}=p_{n}$, and the matrix $\Theta_{{\cal
S}}P=V:=V_{\leq n}({\cal S})$ is our Vandermonde matrix. Therefore, if we can
find a sample set ${\cal S}$ such that the polynomial interpolation problem
with respect to $\Theta_{{\cal S}}$ is poised for $P_{n}$, this would
immediately imply that $|V|\neq 0$. Examples of such situations are provided
in [10, 17], and it turns out that there is a general sufficient condition
which ensures poisedness.

######  Definition 8 ([10]).

A lattice $\mathcal{J}:=\\{x_{1},...,x_{N}\\}$ of $N=p_{n}$ distinct nodes of
$\mathbb{R}^{d}$ is said to satisfy the Geometric Characterization Condition
(GC) (of degree $n$) if corresponding to each node $x_{i}$ there exist $n$
distinct hyperplanes $G_{i1},...,G_{in}$ such that:

  1. 1.

$x_{i}$ does not lie on any of these hyperplanes, and

  2. 2.

all the other nodes in $\mathcal{J}$ lie on at least one of these hyperplanes.

The following result is fundamental for our purposes. We give a short proof
for completeness.

######  Theorem 3 (Theorem 1 in [10]).

Let $\mathcal{J}:=\\{x_{1},...,x_{N}\\}$ be a lattice of $N=p_{n}$ distinct
nodes of $\mathbb{R}^{d}$. If $\mathcal{J}$ satisfies the GC condition as in
Definition 8, then the polynomial interpolation problem with respect to
$\mathcal{J}$ is poised for $P_{n}$.

###### Proof.

Let $Y=\\{y_{i}\\}_{i=1,...,N}$ be an arbitrary vector of values. It can be
directly verified that the following polynomial of total degree $\leq n$
interpolates $Y$ on $\mathcal{J}$:

$P(x)=\sum_{i=1}^{N}y_{i}\prod_{j=1}^{n}{{G_{ij}(x)}\over{G_{ij}(x_{i})}}.$

Since $Y$ was arbitrary, uniqueness follows. ∎

Next, we recall the definition of $n$-th order principal lattice [31].

######  Definition 9.

Let $I_{n}=\\{0,{1\over n},{2\over n},...,1\\}$. Further, denote by
$X_{i}:=e_{i}$ the simplex vertices with $e_{i}$ being the $i$-th element in
standard basis of $\mathbb{R}^{d}$, and furthermore put
$X_{0}:=\mathbf{0}\in\mathbb{R}^{d}$. Then the $n$-th order principal lattice
(corresponding to the coordinate simplex) is defined as

$B(n,d):=\biggl{\\{}x\in\mathbb{R}^{d}\Big{|}x=\sum_{i=0}^{d}\gamma_{i}X_{i},\;\sum_{i=0}^{d}\gamma_{i}=1,\;\gamma_{i}\in
I_{n}\biggr{\\}}.$

Note that $B(n,d)$ contains precisely $p_{n}={n+d\choose n}$ points.

It was proven in [10] Theorem 4, that $B(n,d)$ satisfies the GC condition
which, by Theorem 3, guarantees poisedness [17]. This is almost what we want,
and indeed as we will show, the uniform grid contains a scaled principal
lattice.

######  Lemma 4.

Let $J_{n}:=\\{0,1,...,n\\}$, and put

$A(n,d):=n\cdot
B(n,d)=\biggl{\\{}x\in\mathbb{R}^{d}\Big{|}x=\sum_{i=0}^{d}\gamma_{i}X_{i},\sum_{i=0}^{d}\gamma_{i}=n,\gamma_{i}\in
J_{n}\biggr{\\}}\subset{\cal G}_{n}.$

Then the set $A(n,d)$ satisfies the GC condition of degree $n$.

###### Proof.

The proof is very similar to the one in [10] for $B(n,d)$, and we provide it
for completeness.

Let $\\{\gamma_{k}:\mathbb{R}^{d}\to\mathbb{R}\\}_{k=0,...,d}$ be the
barycentric coordinate functions associated with the simplex
$X_{0},...,X_{d}$. Define the following hyperplanes
$H_{km}:=\\{x|\gamma_{k}(x)=m\\}$, $k=0,...,d$, $m=0,...,n$. These hyperplanes
contain all the nodes in $A(n,d)$.

Now let $x\in A(n,d)$ be a lattice node, and let $\beta_{k}:=\gamma_{k}(x)$.
For each $k=0,..,d$ with $\beta_{k}>0$, the union of the hyperplanes
$\\{H_{kj}\\}_{k=0,...,d}^{j=0,...,\beta_{k}-1}$ doesn’t contain the node $x$.
Since $\sum_{k=0}^{d}\beta_{k}=n$ we have exactly $n$ hyperplanes in the set
$\\{H_{kj}\\}_{k=0,...,d}^{j=0,...,\beta_{k}-1}$.

We still have to show why the union of
$\\{H_{kj}\\}_{k=0,...,d}^{j=0,...,\beta_{k}-1}$ contains all the other nodes.
Let $y\neq x\in A(n,d)$. Again, we have $\sum_{k=0}^{d}\gamma_{k}(y)=n$,
therefore there is at least one $k\in\\{0,...,d\\}$ such that
$\gamma_{k}(y)<\gamma_{k}(x)=\beta_{k}$. ∎

To clarify the proof, we give an example for $d=2,n=2$. For the node $x=(0,1)$
we have $\gamma_{0}(x)=1,\gamma_{1}(x)=0,\gamma_{2}(x)=1$ because
$x=1\cdot(0,0)+0\cdot(1,0)+1\cdot(0,1)$. Following the proof, $H_{00}$ and
$H_{20}$ are the hyperplanes which contain all the nodes except for $x$ (see
figure 1).

Figure 1: Demonstration of the proof of Lemma 4.

######  Remark 5.

For the scaled principal lattice $A(n,d)$, the interpolation polynomial $P(x)$
of degree $\leq n$ for any $\\{y_{a}\\}_{a\in A(n,d)}$ is given by:

$P(x)=\sum_{a\in
A(n,d)}y_{a}\prod_{k=0}^{d}\prod_{m=0}^{\gamma_{k}(a)-1}{{\gamma_{k}(x)-m}\over{\gamma_{k}(a)-m}}.$

We are finally in a position to state and prove the main result of this
section.

######  Theorem 4.

Let ${\cal X}\subset\mathbb{R}^{d}$ be a finite set of $n$ distinct points in
$\mathbb{R}^{d}$, and let $m=\mu({\cal X})$ be its discrete moment order. Let
$\ell$ be the unique integer such that $p_{\ell-1}<n\leq p_{\ell}$. Then for
every $N\geq\ell$ we have

$\operatorname{rank}(V_{\leq m}({\cal G}_{N}))\geq n,$

so the conclusion of Lemma 3 holds for ${\cal S}={\cal G}_{N}$.

###### Proof.

Using Lemma 4 and Theorems 2 and 3, we get

$\displaystyle\big{|}V_{\leq\ell}(A(\ell,d))\big{|}\neq 0$
$\displaystyle\implies\operatorname{rank}(V_{\leq\ell}(A(\ell,d)))={{\ell+d}\choose
d}=p_{\ell}\geq n$ $\displaystyle A(\ell,d)\subset{\cal G}_{\ell}$
$\displaystyle\implies\operatorname{rank}(V_{\leq\ell}({\cal
G}_{\ell}))\geq{{\ell+d}\choose d}\geq n.$

Now notice that since $\operatorname{rank}V_{\leq m}({\cal X})=n$ we must have
$p_{m}\geq n$ and therefore $m\geq\ell$. This implies

$\operatorname{rank}(V_{\leq m}({\cal G}_{N}))\geq\operatorname{rank}(V_{\leq
m}({\cal G}_{\ell})\geq\operatorname{rank}(V_{\leq\ell}({\cal G}_{\ell}))\geq
n,$

finishing the proof. ∎

##  5 Vandermonde matrices in the super-resolution regime

In this section we come back to the super-resolution problem, and consider the
Vandermonde matrices $U(x_{1},\dots,x_{n};{\cal G}_{N})$ in (3) in the single
cluster setting. We prove the multidimensional analogue of [7, Theorem 2.3],
with the caveat that the geometry of the nodes stays fixed as $\text{SRF}\gg
1$. We start by defining the cluster geometry. Naturally, since the sampling
set is an integer grid, we must restrict the nodes to
$\mathbb{T}^{d}:=[-\pi,\pi)^{d}\equiv\left(\mathbb{R}\mod 2\pi\right)^{d}$ to
avoid aliasing.

######  Definition 10.

For $x,y\in\mathbb{T}^{d}$, we denote the wrap-around distance by

$\|x-y\|_{\mathbb{T}^{d}}:=\min_{r\in(2\pi\mathbb{Z})^{d}}\|x-y+r\|_{\infty}.$

######  Definition 11.

Let $\mathcal{X}=\\{x_{1},...,x_{n}\\}\subset\mathbb{T}^{d}$. If for some
$\tau>1$ and $0<\Delta<\pi/\tau$ we have

$\forall x,y\in\mathcal{X},x\neq
y:\quad\Delta\leq\|x-y\|_{\mathbb{T}^{d}}\leq\tau\Delta,$

then $\mathcal{X}$ is said to form an $(\Delta,\tau,n)$-cluster.

In what follows we fix a node set ${\cal
Y}=\\{y_{1},\dots,y_{n}\\}\subset\mathbb{T}^{d}$ such that the (non-wrapped)
minimal distance between any two nodes is

$\rho=\rho({\cal Y}):=\min_{i\neq j}\|y_{i}-y_{j}\|_{\infty}.$

Now let ${\cal X}=\Delta{\cal Y}$ with $\Delta<\frac{1}{2\pi}$, so that ${\cal
X}$ forms an $(\Delta^{\prime},\tau^{\prime},n)$-cluster for
$\Delta^{\prime}=\rho\Delta$ and some $\tau^{\prime}\leq\frac{2\pi}{\rho}$.
(Otherwise, there exist $y,y^{\prime}\in{\cal Y}$ such that
$\|y-y^{\prime}\|_{\mathbb{T}^{d}}\geq\frac{2\pi}{\rho}\Delta^{\prime}=2\pi\Delta$,
a contradiction).

Further, let ${\cal G}_{N}=\\{0,\pm 1,...,\pm N\\}^{d}$ be the uniform
symmetric grid of side length $2N+1$.

Let $U:=U(x_{1},\dots,x_{n};{\cal G}_{N})$. Denote $\epsilon:=N\Delta$. We
have

$U^{*}U=D({\cal X},{\cal G}_{N})=D\Bigl{(}\epsilon\frac{{\cal Y}}{N},{\cal
G}_{N}\Bigr{)}=D_{\epsilon}\Bigl{(}\frac{{\cal Y}}{N},{\cal G}_{N}\Bigr{)}.$

We further define the rescaled Dirichlet kernel

$\displaystyle\mathcal{K}(x,y)$
$\displaystyle=\mathcal{K}_{N}(x,y):=\frac{1}{(2N)^{d}}\mathcal{D}_{{\cal
G}_{N}}(x/N,y/N)$ $\displaystyle\implies K_{\epsilon}({\cal Y})$
$\displaystyle=\frac{1}{(2N)^{d}}D_{\epsilon}\Bigl{(}\frac{{\cal Y}}{N},{\cal
G}_{N}\Bigr{)}=D_{N}({\cal X}).$

We have the following result.

######  Theorem 5.

With the above notations, let $m=\mu({\cal Y})=\mu({\cal X})$ be the discrete
moment order of ${\cal Y}$.

For every $\epsilon_{0}>0,\alpha\in(0,1)$ there exists $N_{0}$ such that for
all $N\geq N_{0}$ and $\alpha\epsilon_{0}<N\Delta<\epsilon_{0}$, the
eigenvalues of $D_{N}$ split into $m+1$ groups

$\lambda_{0,0}(D_{N})\asymp{(N\Delta)}^{0},\quad\\{\lambda_{1,j}(D_{N})\\}_{j=1}^{t_{1}}\asymp{(N\Delta)}^{2},...,\\{\lambda_{m,j}(D_{N})\\}_{j=1}^{t_{m}}\asymp{(N\Delta)}^{2m},$

where $t_{k}:=\operatorname{rank}(V_{\leq k}({\cal
X}))-\operatorname{rank}(V_{\leq k-1}({\cal X}))$, $k=1,\dots,m$.

###### Proof.

With the identification $D_{N}({\cal X})=K_{\epsilon}({\cal Y})$ as above, we
would like to apply Theorem 1 and Corollary 1 to the kernel
$\mathcal{K}=\mathcal{K}_{N}$. However, $\mathcal{K}$ depends on $N$, so we
need to make sure that all the estimates in the proofs are uniform in $N$.

######  Definition 12 (sinc kernel).

Let $\operatorname{sinc}(x):=\begin{cases}\frac{\sin(x)}{x},&\text{if }x\neq
0,\\\ 1,&\text{if }x=0.\end{cases}$. Given ${\cal X}\in\mathbb{R}^{d}$, the
corresponding sinc kernel matrix is defined as follows:

$\operatorname{Sinc}({\cal
X}):=\biggl{[}\prod_{i=1}^{d}\operatorname{sinc}(x_{i}-x_{i}^{\prime})\biggr{]}_{x,x^{\prime}\in{\cal
X}}$ (12)

For $\epsilon>0$, the flat limit version of the $\operatorname{sinc}$ kernel
matrix is:

$\operatorname{Sinc}_{\epsilon}({\cal
X}):=\biggl{[}\prod_{i=1}^{d}\operatorname{sinc}(\epsilon x_{i}-\epsilon
x_{i}^{\prime})\biggr{]}_{x,x^{\prime}\in{\cal X}}$ (13)

Outline of the proof: we aim to prove that the eigenvalues of $D_{N}$ converge
to $\operatorname{Sinc}$ kernel when $N\to\infty$. The scaling of the
eigenvalues of $\operatorname{Sinc}$ can be obtained using Corollary 1 (the
constants $\tilde{\lambda}_{k,j}$ don’t depend on $N$). First we show element-
wise convergence of $D_{N}$ to $\operatorname{sinc}$ and equicontinuity of
$D_{N}$ in $\epsilon<\epsilon_{0}$ to show uniform convergence of $D_{N}$ to
$\operatorname{sinc}$ (uniform in $\epsilon<\epsilon_{0}$ and independent of
$N$). Then using the Bauer–Fike Theorem, we prove convergence of the
eigenvalues.

We start with some definitions.

######  Theorem 6 (The Bauer–Fike Theorem).

For a diagonalizable matrix $A$ let $E$ be it’s eigenvectors matrix,
$A=E\Lambda{E}^{-1}$ where $\Lambda$ is a diagonal matrix. Let $\mu$ be an
eigenvalue of $A_{N}$, then there exists $\lambda\in\Lambda(A)$ such that

$|\lambda-\mu|\leq\kappa_{p}(E)\|A-A_{N}\|_{p}$ (14)

where $\kappa_{p}(E):=\|E\|_{p}\|{E}^{-1}\|_{p}$ is the condition number of
$E$ associated with the norm $\|\cdot\|_{p}$.

######  Definition 13 (Element-wise convergence).

Let $\\{A_{N}\\}_{N=1}^{\infty}$, where $A_{N}\in\mathbb{R}^{s\times s}$, be a
sequence of matrices, and let $A\in\mathbb{R}^{s\times s}$ be a matrix. The
sequence $\\{A_{N}\\}_{N=1}^{\infty}$ is said to converge element-wise to $A$
if the following condition is satisfied:

For every $\delta>0$, there exists a positive integer $M(\delta)$ such that
for all $N\geq M(\delta)$, the inequality

$\big{|}(A_{N})_{i,j}-(A)_{i,j}\big{|}\leq\delta$

holds for all $1\leq i,j\leq s$.

######  Definition 14 (Equicontinuity).

Let $\\{f_{n}(x)\\}$ be a sequence of functions. $\\{f_{n}(x)\\}$ is
equicontinuous if for every $\delta>0$, there exists $\nu>0$ (independent of
$n$) such that for all $n$ and all $x,y\in[a,b]$ with $|x-y|<\nu$ we have
$|f_{n}(x)-f_{n}(y)|<\delta$.

######  Definition 15 (Uniform convergence).

A sequence of functions $\\{f_{n}(x)\\}$ converges uniformly to $f(x)$ on $S$
if for every $\epsilon>0$, there exists an integer $N$ such that for all
$n\geq N$ and for all $x\in S$, $|f_{n}(x)-f(x)|<\epsilon$.

######  Theorem 7.

Let $C(J)$ be the space of continuous functions on a closed interval
$J\subset\mathbb{R}$. A sequence in $C(J)$ is uniformly convergent if and only
if it is equicontinuous and converges pointwise to a function.

###### Proof.

The first direction is an immediate corollary of Theorem 7.25 in [35]. Since
the sequence converges pointwise, we can replace the subsequence in the proof
of (b) by the sequence itself. The second direction is Theorem 7.24. ∎

######  Lemma 5.

For every $\epsilon>0$ and ${\cal X}\in\mathbb{R}^{d}$, the Dirichlet kernel
matrix $D_{N}({\cal X})$ converges element-wise to the $\operatorname{sinc}$
kernel matrix $\operatorname{Sinc}_{\epsilon}(\mathcal{Y})$ when $N\to\infty$.

###### Proof.

$\big{(}D_{N}({\cal
X})\big{)}_{i,j}=\frac{1}{(2N)^{d}}\big{(}D_{\epsilon}\Bigl{(}\frac{{\cal
Y}}{N},{\cal G}_{N}\Bigr{)}\big{)}_{i,j}=\frac{1}{(2N)^{d}}\mathcal{D}_{{\cal
G}_{N}}(\epsilon({y_{i}\over N}-{y_{j}\over
N}))=\frac{1}{(2N)^{d}}\sum_{\omega\in{\cal
G}_{N}}\exp\left(\imath\langle\epsilon({y_{i}\over N}-{y_{j}\over
N}),\omega\rangle\right)$ (15)

Denote by $y:=y_{i}-y_{j}$ and $y:=(y^{(1)},\dots,y^{(d)})$. First, we
consider the case where $d=1$:

$\displaystyle\frac{1}{2N}\sum_{k=-N}^{N}\exp(\imath\epsilon k{y\over
N})=\frac{1}{2N}(1+2\sum_{k=1}^{N}\cos{(\epsilon k{y\over
N})})=\frac{1}{2N}\frac{\sin((N+{1\over 2}){\epsilon y\over
N})}{\sin({\epsilon y\over 2N})}$
$\displaystyle\to_{N\to\infty}\frac{1}{2N}\frac{\sin((1+{1\over 2N})\epsilon
y)}{{\epsilon y\over 2N}}\to_{N\to\infty}{\sin(\epsilon y)\over\epsilon
y}=\operatorname{sinc}(\epsilon y).$

For $d>1$, using the above, we get

$\displaystyle\frac{1}{(2N)^{d}}\sum_{\omega\in{\cal
G}_{N}}\exp\left(\imath\langle\epsilon({y\over
N}),\omega\rangle\right)=\prod_{i=1}^{d}\frac{1}{2N}\sum_{k=-N}^{N}\exp(\imath\epsilon
k{y^{(i)}\over N})\to_{N\to\infty}\prod_{i=1}^{d}\operatorname{sinc}(\epsilon
y^{(i)}).$

∎

We proved that each entry of $D_{N}({\cal X})$ converges point-wise to each
entry of $\operatorname{Sinc}_{\epsilon}(\mathcal{Y})$. Now we show that this
convergence is uniform for $\epsilon\in[0,\epsilon_{0}]$ for a given
$\epsilon_{0}$. We show that using Theorem 7. First we need to show that the
sequence $\\{\big{(}D_{N}({\cal X})\big{)}_{i,j}\\}$ is equicontinuous.

For fixed $y\in\mathcal{Y}$, let
$f_{N}(\epsilon):=\frac{1}{2N}\sum_{k=-N}^{N}\exp(\imath\epsilon k{y\over N})$
and $f(\epsilon):=\operatorname{sinc}(\epsilon y)$. We prove the following
result.

######  Lemma 6.

Let $\epsilon_{0}>0$ be given. For every $\delta>0$ there exists $\nu>0$
(independent of $N$) such that for all $N$ and all
$\epsilon,\epsilon^{\prime}\in[0,\epsilon_{0}]$ with
$|\epsilon-\epsilon^{\prime}|<\nu$ we have
$|f_{n}(\epsilon)-f_{n}(\epsilon^{\prime})|<\delta$.

###### Proof.

Let $\delta>0$. we want to find $\nu(\delta)>0$ such that for all $N$ and
$\epsilon,\epsilon^{\prime}\in[0,\epsilon_{0}]$ with
$|\epsilon-\epsilon^{\prime}|<\nu$ we have that
$|f_{N}(\epsilon)-f_{N}(\epsilon^{\prime})|<\delta$. Using that
$|1-e^{ix}|=2|\sin({x\over 2})|$ and $\sin(x)\leq x$ for $x\in[0,{\pi\over
2}]$ we get

$\displaystyle|f_{N}(\epsilon)-f_{N}(\epsilon^{\prime})|$
$\displaystyle=|\frac{1}{2N}\sum_{k=-N}^{N}\exp(\imath\epsilon k{y\over
N})-\frac{1}{2N}\sum_{k=-N}^{N}\exp(\imath\epsilon^{\prime}k{y\over N})|$
$\displaystyle=|\frac{1}{2N}\sum_{k=-N}^{N}\exp(\imath\epsilon k{y\over
N})(1-\exp(\imath(\epsilon-\epsilon^{\prime})k{y\over N}))|$
$\displaystyle\leq\frac{1}{2N}\sum_{k=-N}^{N}|1-\exp(\imath(\epsilon-\epsilon^{\prime})k{y\over
N})|$
$\displaystyle=\frac{1}{2N}\sum_{k=-N}^{N}2|\sin((\epsilon-\epsilon^{\prime})k{y\over
2N})|$
$\displaystyle\leq\frac{2}{2N}\sum_{k=0}^{N}2k|\epsilon-\epsilon^{\prime}|{y\over
2N}={1\over N^{2}}{N+N^{2}\over 2}\nu y$ $\displaystyle=({1\over 2N}+{1\over
2})y\nu\leq y\nu.$

Thus we choose $\nu:={\delta\over y}$. The same holds for
$f_{N}(\epsilon):=\frac{1}{(2N)^{d}}\sum_{\omega\in{\cal
G}_{N}}\exp\left(\imath\langle\epsilon({y\over N}),\omega\rangle\right)$. ∎

Using Lemma 5 and 6 together with Theorem 7 we get the following Corollary.

######  Corollary 2.

Let $\epsilon_{0}$ be given. For every $\delta>0$ there exists
$M=M(\epsilon_{0},\delta)$ such that for all $N>M$

$\big{|}\big{(}D_{N}({\cal
X})\big{)}_{i,j}-\big{(}\operatorname{Sinc}_{\epsilon}(\mathcal{Y})\big{)}_{i,j}\big{|}\leq\delta,\quad\forall\epsilon\in[0,\epsilon_{0}].$
(16)

###### Proof.

We take $\nu$ in Definition 14 to be the minimum of all the
$\nu_{i,j}^{(\ell)}$ according to the proof of Lemma 6, with $\ell=1,\dots,d$.
Note that for $y_{i,j}^{(\ell)}=0$ we have $\nu_{i,j}^{(\ell)}=\infty$. ∎

######  Corollary 3 (Uniform convergence of eigenvalues).

Let $\epsilon_{0}>0$ be given. For every $\delta^{\prime}>0$ there exists
$M=M(\epsilon_{0},\delta^{\prime})$ s.t. for all $N>M$ there is an eigenvalue
$\mu_{k}^{(N)}$ of $D_{N}({\cal X})$ with

$|\lambda^{(\epsilon)}_{k}-\mu_{k}^{(N)}|\leq\delta^{\prime},\quad\forall\epsilon\in[0,\epsilon_{0}]$
(17)

where $\lambda^{(\epsilon)}_{k}$ is an eigenvalue of
$\operatorname{Sinc}_{\epsilon}(\mathcal{Y})$.

###### Proof.

Since $\operatorname{Sinc}_{\epsilon}(\mathcal{Y})$ is a symmetric matrix, the
eigenvectors matrix $E$ is orthogonal in the eigenvalue decomposition of
$\operatorname{Sinc}_{\epsilon}(\mathcal{Y})$ i.e
$\operatorname{Sinc}_{\epsilon}(\mathcal{Y})=E\Lambda{E}^{-1}$. Thus we have
$\kappa_{2}(E)=1$. Using the Bauer–Fike theorem 6 with Corollary 2 and
$\delta:={\delta^{\prime}\over n\sqrt{n}}$, we get:

$\displaystyle|\lambda-\mu|\leq\|\operatorname{Sinc}_{\epsilon}(\mathcal{Y})-D_{N}({\cal
X})\|_{2}$
$\displaystyle\leq\sqrt{n}\|\operatorname{Sinc}_{\epsilon}(\mathcal{Y})-D_{N}({\cal
X})\|_{\infty}$
$\displaystyle=\sqrt{n}\max_{i}\sum_{j=1}^{n}\big{|}(\operatorname{Sinc}_{\epsilon}(\mathcal{Y}))_{i,j}-(D_{N}({\cal
X}))_{i,j}\big{|}$ $\displaystyle\leq\sqrt{n}n\delta=\delta^{\prime}.$

∎

To complete the proof of Theorem 5, we choose
$\delta^{\prime}=O\biggl{(}(\alpha\epsilon_{0})^{2m}\biggr{)}$ in (17) and
apply Corollary 1. The only remaining question is to show that $\det(RWR)\neq
0$. We will show that the limit $N\to\infty$ the rank of
$W=\lim_{N\to\infty}\frac{1}{(2N)^{d}}W_{N}$ is full, and as a result,
$\det\left(\tilde{R}W\tilde{R}^{T}\right)$ is uniformly bounded away from
zero. By definition of the Wronskian, we write

$\displaystyle\frac{1}{(2N)^{d}}W_{N}$
$\displaystyle=\frac{1}{(2N)^{d}}\bigg{[}{{\mathcal{D}_{{\cal
G}_{N}}^{(\alpha,\beta)}({x\over N},{y\over
N})|_{x,y=0}}\over{\alpha!\beta!}}\bigg{]}_{\alpha,\beta\in\mathbb{P}_{m}}=\bigg{[}\frac{1}{(2N)^{d}}\sum_{\omega\in{\cal
G}_{N}}{1\over
N^{|\alpha|+|\beta|}}{{(i\omega)^{\alpha}(-i\omega)^{\beta}}\over{\alpha!\beta!}}\bigg{]}_{\alpha,\beta\in\mathbb{P}_{m}}$
$\displaystyle=\bigg{[}\frac{1}{(2N)^{d}}\sum_{\omega=(\omega_{1},\dots,\omega_{d})\in{\cal
G}_{N}}{{(i)^{|\alpha|}(-i)^{|\beta|}}\over{\alpha!\beta!}}\big{(}{\omega_{1}\over
N}\big{)}^{\alpha_{1}+\beta_{1}}\cdot...\cdot\big{(}{\omega_{d}\over
N}\big{)}^{\alpha_{d}+\beta_{d}}\bigg{]}_{\alpha,\beta\in\mathbb{P}_{m}}$
$\displaystyle=\bigg{[}\prod_{j=1}^{d}\sum_{k=-N}^{N}\frac{1}{2N}\frac{(i)^{\alpha_{j}}(-i)^{\beta_{j}}}{\alpha_{j}!\beta_{j}!}\big{(}{k\over
N}\big{)}^{\alpha_{j}+\beta_{j}}\bigg{]}_{\alpha,\beta\in\mathbb{P}_{m}}.$

We get that each entry of $\frac{1}{(2N)^{d}}W_{N}$ is a Riemann sum, thus as
$N\to\infty$ we get a multiple integral:

$\lim_{N\to\infty}\prod_{j=1}^{d}\sum_{k=-N}^{N}\frac{1}{2N}\frac{(i)^{\alpha_{j}}(-i)^{\beta_{j}}}{\alpha_{j}!\beta_{j}!}\biggl{(}{k\over
N}\biggr{)}^{\alpha_{j}+\beta_{j}}={1\over
2^{d}}\int_{[-1,1]^{d}}{{(i)^{|\alpha|}(-i)^{|\beta|}}\over{\alpha!\beta!}}x^{\alpha+\beta}\,dx$

Thus we get:

$W=\lim_{N\to\infty}\frac{1}{(2N)^{d}}W_{N}=\bigg{[}{1\over
2^{d}}\int_{[-1,1]^{d}}{{(i)^{|\alpha|}(-i)^{|\beta|}}\over{\alpha!\beta!}}x^{\alpha+\beta}\,dx\bigg{]}_{\alpha,\beta\in\mathbb{P}_{m}}.$

It is immediately seen that $W$ is none other than the Gram matrix for the
monomials
$\\{\boldsymbol{x}^{\boldsymbol{\alpha}}\\}_{\boldsymbol{\alpha}\in\mathbb{P}_{m}}$:

$W_{\alpha,\beta}=\langle x^{\alpha},x^{\beta}\rangle={1\over
2^{d}}\int_{[-1,1]^{d}}{{(i)^{|\alpha|}(-i)^{|\beta|}}\over{\alpha!\beta!}}x^{\alpha+\beta}\,dx.$

Since the set of monomials is linearly independent, we have that $W$ is
positive definite. Thus $W$ is full rank, and as a result,
$\det\left(\tilde{R}W\tilde{R}^{T}\right)>0$. This completes the proof of
Theorem 5. ∎

######  Remark 6.

During the final stages of preparation of our manuscript, we became aware that
Weilin Li is working on a related topic [24]. Relative to our results, more
explicit conditions and constants are derived for the scaling of the smallest
singular value of the multivariate Vandermonde matrix (we consider the entire
spectrum); multi-cluster geometry is investigated; only cube and spherical
sampling sets are considered; the decay rate in the case of uniform dilations
of generic sets is not sharp; finally, a form of geometric characterization
condition is used to describe the node geometry (as opposed to the sequence
$\\{t_{k}\\}$ in our work).

##  6 Numerical Experiments

In this section, we describe some numerical example to validate our results
and investigate their applicability to the problem of super-resolution.

###  6.1 Dirichlet kernel eigenvalue scaling

(a) Parabola geometry.

(b) General position geometry.

Figure 2: Eigenvalues of Dirichlet kernel with samples on the grid for the
different node geometries. The $x$-axis is $\epsilon$ denoting the smallest
distance between the nodes.

We start with the numerical study of the eigenvalues of the Dirichlet matrix
$D_{\epsilon}({\cal X},{\cal G}_{n})$, confirming the scaling predicted by
Lemma 3. Let $d=2$, $n=6$ and $\mathcal{X}=\\{(x_{i},y_{i})\\}_{i=1}^{n}$. We
checked two scenarios:

  1. 1.

All points satisfy $y_{i}=x_{i}^{2}$. We get that $\operatorname{rank}(V_{\leq
3})=6$, where

$\displaystyle V_{\leq 3}$
$\displaystyle=\left[\begin{array}[]{c|cc|ccc|cccc}1&x_{1}&y_{1}&x_{1}^{2}&x_{1}y_{1}&y_{1}^{2}&x_{1}^{3}&x_{1}^{2}y_{1}&x_{1}y_{1}^{2}&y_{1}^{3}\\\
1&x_{2}&y_{2}&x_{2}^{2}&x_{2}y_{2}&y_{2}^{2}&x_{2}^{3}&x_{2}^{2}y_{2}&x_{2}y_{2}^{2}&y_{2}^{3}\\\
1&x_{3}&y_{3}&x_{3}^{2}&x_{3}y_{3}&y_{3}^{2}&x_{3}^{3}&x_{3}^{2}y_{3}&x_{3}y_{3}^{2}&y_{3}^{3}\\\
1&x_{4}&y_{4}&x_{4}^{2}&x_{4}y_{4}&y_{4}^{2}&x_{4}^{3}&x_{4}^{2}y_{4}&x_{4}y_{4}^{2}&y_{4}^{3}\\\
1&x_{5}&y_{5}&x_{5}^{2}&x_{5}y_{5}&y_{5}^{2}&x_{5}^{3}&x_{5}^{2}y_{5}&x_{5}y_{5}^{2}&y_{5}^{3}\\\
1&x_{6}&y_{6}&x_{6}^{2}&x_{6}y_{6}&y_{6}^{2}&x_{6}^{3}&x_{5}^{2}y_{5}&x_{6}y_{6}^{2}&y_{6}^{3}\end{array}\right]$
$\displaystyle=\left[\begin{array}[]{c|cc|ccc|cccc}1&x_{1}&x^{2}_{1}&x_{1}^{2}&x_{1}^{3}&x_{1}^{4}&x_{1}^{3}&x_{1}^{3}&x_{1}^{4}&x_{1}^{6}\\\
1&x_{2}&x^{2}_{2}&x_{2}^{2}&x_{2}^{3}&x_{2}^{4}&x_{2}^{3}&x_{2}^{3}&x_{2}^{4}&x_{2}^{6}\\\
1&x_{3}&x^{2}_{3}&x_{3}^{2}&x_{3}^{3}&x_{3}^{4}&x_{3}^{3}&x_{3}^{3}&x_{3}^{4}&x_{3}^{6}\\\
1&x_{4}&x^{2}_{4}&x_{4}^{2}&x_{4}^{3}&x_{4}^{4}&x_{4}^{3}&x_{4}^{3}&x_{4}^{4}&x_{4}^{6}\\\
1&x_{5}&x^{2}_{5}&x_{5}^{2}&x_{5}^{3}&x_{5}^{4}&x_{5}^{3}&x_{5}^{3}&x_{5}^{4}&x_{5}^{6}\\\
1&x_{6}&x^{2}_{6}&x_{6}^{2}&x_{6}^{3}&x_{6}^{4}&x_{6}^{3}&x_{5}^{3}&x_{6}^{4}&x_{6}^{6}\end{array}\right]=\begin{bmatrix}V_{0}&|&V_{1}&|&V_{2}&|&V_{3}\end{bmatrix}$

We get that

$\displaystyle t_{0}$ $\displaystyle=\operatorname{rank}(V_{\leq 0})=1,$
$\displaystyle t_{1}$ $\displaystyle=\operatorname{rank}(V_{\leq
1})-rank(V_{\leq 0})=3-1=2,$ $\displaystyle t_{2}$
$\displaystyle=\operatorname{rank}(V_{\leq 2})-rank(V_{\leq 1})=5-3=2,$
$\displaystyle t_{3}$ $\displaystyle=\operatorname{rank}(V_{\leq
3})-\operatorname{rank}(V_{\leq 2})=6-5=1.$

The results in Figure 2(a) confirm the scaling of the eigenvalues.

  2. 2.

The set of nodes in ${\cal X}$ are in general positions. We have
$\operatorname{rank}(V_{\leq 2})=6$, where

$\displaystyle V_{\leq
2}=\left[\begin{array}[]{c|cc|ccc}1&x_{1}&y_{1}&x_{1}^{2}&x_{1}y_{1}&y_{1}^{2}\\\
1&x_{2}&y_{2}&x_{2}^{2}&x_{2}y_{2}&y_{2}^{2}\\\
1&x_{3}&y_{3}&x_{3}^{2}&x_{3}y_{3}&y_{3}^{2}\\\
1&x_{4}&y_{4}&x_{4}^{2}&x_{4}y_{4}&y_{4}^{2}\\\
1&x_{5}&y_{5}&x_{5}^{2}&x_{5}y_{5}&y_{5}^{2}\\\
1&x_{6}&y_{6}&x_{6}^{2}&x_{6}y_{6}&y_{6}^{2}\end{array}\right]=\begin{bmatrix}V_{0}&|&V_{1}&|&V_{2}\end{bmatrix}.$

We get that

$\displaystyle t_{0}$ $\displaystyle=\operatorname{rank}(V_{\leq 0})=1,$
$\displaystyle t_{1}$ $\displaystyle=\operatorname{rank}(V_{\leq
1})-\operatorname{rank}(V_{\leq 0})=3-1=2,$ $\displaystyle t_{2}$
$\displaystyle=\operatorname{rank}(V_{\leq 2})-\operatorname{rank}(V_{\leq
1})=6-3=3.$

The results in Figure 2(b) confirm the scaling of the eigenvalues also in this
case.

###  6.2 Super-resolution

In this section, we consider the super-resolution problem for multidimensional
sparse measures described in the Introduction, as it was the original
motivation for our investigations. We reconstruct measures of the form (1)
from the noisy measurements (2). We consider the asymptotic regime of constant
$N$ and $\Delta\to 0$.

####  6.2.1 Local Stability

Analogous to some previous studies in the one-dimensional case, e.g. [2, 4],
we consider the “local stability” of the problem to be well-approximated by
the error incurred by the nonlinear least squares (NLS) method applied to the
noisy measurements:

$\left\\{\hat{x}_{j}^{NLS},\hat{\alpha}_{j}^{NLS}\right\\}_{j=1}^{n}=\arg\min_{\hat{\alpha},\hat{x}}\frac{1}{2}\sum_{k\in{\cal
G}_{N}}\left|\hat{f}(k)-\sum_{j=1}^{n}\hat{\alpha}_{j}e^{\imath\left\langle
k,\hat{x}_{j}\right\rangle}\right|^{2}.$

The initial values for the NLS method are taken to be the true parameter
values (which are clearly unknown in practice). Our implementation of the NLS
method is based on the Levenberg-Marquardt algorithm, with the complex
residuals converted to real residuals by concatenating the real and imaginary
parts. The results are shown in Figures 3 and 4. To measure the accuracy, we
fitted the scaled error

$\kappa:=\max_{j=1,\dots,n}\|\hat{x}_{j}^{NLS}-x_{j}\|_{\infty}/\sigma$

to the curve $\kappa\sim\Delta^{m}$, where $\Delta$ is the minimal distance
between the nodes. It can be seen that the accuracy crucially depends on the
geometry, and is indeed best for the random configuration, and worst for the
single-line configuration (for which the scaling is precisely
$\kappa\sim\Delta^{2n-2}$ for any $d$, as established in previous works on the
subject).

Figure 3: NLS reconstruction, general position/line/parabola, $d=2$.

Figure 4: NLS reconstruction, general position/line/parabola, $d=3$.

####  6.2.2 2D-ESPRIT

We have implemented the 2D ESPRIT algorithm according to [36, Section IV-A].
The experimental setup was as follows:

  1. 1.

The number of nodes varied from $n=2$ to $n=6$.

  2. 2.

The geometry was either random, a line $y=5x$, or a parabola $y=10x^{2}$.

  3. 3.

The ESPRIT algorithm parameters were set as follows:
$M_{1}=M_{2}=40,\;L_{1}=L_{2}=10,\;\beta_{1}=\beta_{2}=0.5$.

  4. 4.

The noise level was set to $\sigma=10^{-20}$.

The results are shown in Figure 5. As for the NLS method, we fitted the scaled
error

$\kappa:=\max_{j=1,\dots,n}\|\hat{x}_{j}^{ESPRIT}-x_{j}\|_{\infty}/\sigma$

to the curve $\kappa\sim\Delta^{m}$, where $\Delta$ is the minimal distance
between the nodes. As in the NLS case, it can be seen that the accuracy
crucially depends on the geometry, providing best scaling for the random
configuration, and worst scaling for the single-line configuration.

Figure 5: Accuracy of 2D ESPRIT for different geometries.

##  Appendix A Taylor Expansion

Let $\mathcal{K}\in\mathcal{C}^{(r,r)}(\Omega)$ where $\Omega$ is an open
neighborhood of $0$ and $r\geq m+1$, where $m=\mu({\cal X})$ is the discrete
moment order of ${\cal X}$ (Definition 5). We recall the Maclaurin expansion
for bivariate functions from [1, section 5]. For $\epsilon x,\epsilon y$ such
that $[0,\epsilon x],[0,\epsilon y]\subset\Omega$ (where $[0,x]$ is a line
segment from $0$ to $x$), we have:

$\displaystyle\mathcal{K}_{\epsilon}(x,y)=\mathcal{K}(\epsilon x,\epsilon y)$
$\displaystyle=\sum_{\alpha,\beta\in\mathbb{P}_{m}}\frac{{(\epsilon
x)}^{\alpha}{(\epsilon
y)}^{\beta}}{\alpha!\beta!}\mathcal{K}^{(\alpha,\beta)}(0,0)$
$\displaystyle+\sum_{\alpha\in\mathbb{P}_{m},\beta\in\mathbb{H}_{m+1}}\frac{{(\epsilon
x)}^{\alpha}{(\epsilon
y)}^{\beta}}{\alpha!\beta!}\mathcal{K}^{(\alpha,\beta)}(0,\theta_{\epsilon
y,\alpha}\epsilon
y)+\sum_{\alpha\in\mathbb{H}_{m+1},\beta\in\mathbb{P}_{m}}\frac{{(\epsilon
x)}^{\alpha}{(\epsilon
y)}^{\beta}}{\alpha!\beta!}\mathcal{K}^{(\alpha,\beta)}(\eta_{\epsilon
x,\beta}\epsilon x,0)$
$\displaystyle+\sum_{\alpha,\beta\in\mathbb{H}_{m+1}}\frac{{(\epsilon
x)}^{\alpha}{(\epsilon
y)}^{\beta}}{\alpha!\beta!}\mathcal{K}^{(\alpha,\beta)}(\psi_{\epsilon
x,\epsilon y}\epsilon x,\xi_{\epsilon x,\epsilon y}\epsilon y),$

where $\\{\theta_{\epsilon y,\alpha}\\}_{\alpha\in\mathbb{P}_{m}}\subset[0,1]$
depend on $\epsilon y$, $\\{\eta_{\epsilon
x,\beta}\\}_{\beta\in\mathbb{P}_{m}}\subset[0,1]$ depend on $\epsilon x$ and
$\psi_{\epsilon x,\epsilon y},\xi_{\epsilon x,\epsilon y}$ depend on both
$\epsilon x$ and $\epsilon y$.

Let $\mathbb{P}_{m}=\\{\alpha_{1},\dots,\alpha_{p_{m}}\\}$. Then we can write
the above expansion as follows:

$\displaystyle\mathcal{K}(\epsilon x,\epsilon y)=\big{[}{(\epsilon
x)}^{\alpha_{1}},\dots,{(\epsilon x)}^{\alpha_{p_{m}}}\big{]}W_{\leq
m}{\big{[}{(\epsilon y)}^{\alpha_{1}},\dots,{(\epsilon
y)}^{\alpha_{p_{m}}}\big{]}}^{T}$
$\displaystyle+\epsilon^{m+1}\big{[}{(\epsilon
x)}^{\alpha_{1}},\dots,{(\epsilon
x)}^{\alpha_{p_{m}}}\big{]}\mathbf{w}_{1,y}(\epsilon)+\epsilon^{m+1}{\mathbf{w}_{2,x}(\epsilon)}^{T}{\big{[}{(\epsilon
y)}^{\alpha_{1}},\dots,{(\epsilon y)}^{\alpha_{p_{m}}}\big{]}}^{T}$
$\displaystyle+\epsilon^{2(m+1)}w_{3,x,y},$
$\displaystyle\mathbf{w}_{1,y}(\epsilon)={\begin{bmatrix}\sum_{\beta\in\mathbb{H}_{m+1}}\frac{\mathcal{K}^{(\alpha_{1},\beta)}(0,\theta_{\epsilon
y,\alpha_{1}}\epsilon
y)}{\alpha_{1}!\beta!}y^{\beta},\dots,\sum_{\beta\in\mathbb{H}_{m+1}}\frac{\mathcal{K}^{(\alpha_{p_{m}},\beta)}(0,\theta_{\epsilon
y,\alpha_{p_{m}}}\epsilon
y)}{\alpha_{p_{m}}!\beta!}y^{\beta}\end{bmatrix}}^{T}$
$\displaystyle\mathbf{w}_{2,x}(\epsilon)=\begin{bmatrix}\sum_{\alpha\in\mathbb{H}_{m+1}}\frac{\mathcal{K}^{(\alpha,\beta)}(\eta_{\epsilon
x,\beta}(\epsilon
x),0)}{\alpha!\beta!}x^{\alpha}\end{bmatrix}_{\beta\in\mathbb{P}_{m}}^{T}$
$\displaystyle
w_{3,x,y}(\epsilon)=\sum_{\alpha,\beta\in\mathbb{H}_{m+1}}\frac{\mathcal{K}^{(\alpha,\beta)}(\psi_{\epsilon
x,\epsilon y}(\epsilon x),\xi_{\epsilon x,\epsilon y}(\epsilon
y))}{\alpha!\beta!}{(\epsilon x)}^{\alpha}{(\epsilon y)}^{\beta}$

where
$\mathbf{w}_{2,x},\mathbf{w}_{1,y}:[0,\epsilon_{0}]\rightarrow\mathbb{R}^{p_{m}}$
are bounded vector functions with $\epsilon_{0}>0$ such that
$\epsilon_{0}x,\epsilon_{0}y\in\Omega$, $\epsilon\in[0,\epsilon_{0}]$ and
$w_{3,x,y}$ is bounded.

Let $\epsilon_{0}>0$ such that $\epsilon_{0}x_{i}\in\Omega$ for $x_{i}\in{\cal
X}$. Thus for $0\leq\epsilon\leq\epsilon_{0}$, the scaled kernel matrix
$K_{\epsilon}$ admits the following expansion:

$K_{\epsilon}=V_{\leq m}\Delta_{m}W\Delta_{m}V_{\leq
m}^{T}+\epsilon^{m+1}(V_{\leq
m}\Delta_{m}W_{1}(\epsilon)+W_{2}(\epsilon)\Delta_{m}V_{\leq
m}^{T})+\epsilon^{2{(m+1)}}W_{3}(\epsilon)$ (18)

where
$W_{1}(\epsilon):=\big{[}\mathbf{w}_{1,x_{1}}(\epsilon),\dots,\mathbf{w}_{1,x_{n}}(\epsilon)\big{]}$,
$W_{2}(\epsilon):={\big{[}\mathbf{w}_{2,x_{1}}(\epsilon),\dots,\mathbf{w}_{2,x_{n}}(\epsilon)\big{]}}^{T}$
and $W_{3}(\epsilon)=\big{[}w_{3,x_{i},x_{j}}(\epsilon)\big{]}_{i,j=1}^{n,n}$.

##  Appendix B Auxiliary Lemmas

######  Lemma 7.

Let $D\in\mathbb{R}^{k\times k}$ be a $k\times k$ diagonal matrix and $P$ a
permutation matrix corresponding to a permutation $\sigma$, then

$PDP^{T}=\operatorname{diag}\\{d_{\sigma(1)},\dots,d_{\sigma(k)}\\},\quad
D=\operatorname{diag}\\{d_{1},\dots,d_{k}\\}.$

###### Proof.

If $P$ is a permutation matrix corresponding to a permutation $\sigma$, then
for every vector $v\in\mathbb{R}^{k}$

$Pv=\begin{bmatrix}v_{\sigma(1)}\\\ \vdots\\\ v_{\sigma(k)}\end{bmatrix}$

and then for every matrix $A=\begin{bmatrix}R_{1}\\\ \vdots\\\
R_{k}\end{bmatrix}\in\mathbb{R}^{k\times k}$, we have

$PA=\begin{bmatrix}R_{\sigma(1)}\\\ \vdots\\\ R_{\sigma(k)}\end{bmatrix}.$

Now let $A=D=\operatorname{diag}\\{d_{1},\dots,d_{k}\\}$, and take
$j=1,\dots,k$. Then $R_{\sigma(j)}=d_{\sigma(j)}e_{\sigma(j)}^{T}$. Therefore

$R_{\sigma(j)}P^{T}=d_{\sigma(j)}\left(Pe_{\sigma(j)}\right)^{T}=d_{\sigma(j)}e_{j}^{T}.$

Stacking all the rows we get
$PDP^{T}=\operatorname{diag}\\{d_{\sigma(1)},\dots,d_{\sigma(k)}\\}$. ∎

######  Lemma 8.

Let $M,B\in\mathbb{R}^{n\times n}$. For $\epsilon>0$ we have

$\det(M+\epsilon B)=\det(M)+O(\epsilon),\quad\epsilon\ll 1.$ (19)

###### Proof.

Recall formula (0.8.12.3) for sum of matrices in [19]:

$\det(M+\epsilon
B)=\sum_{k=0}^{n}\epsilon^{k}\operatorname{tr}(\operatorname{adj}_{k}(M)C_{k}(B)),$
(20)

where $\operatorname{tr}(\cdot)$ is the trace of matrix,
$\operatorname{adj}_{k}(\cdot)$ is the $k$-th adjugate of matrix (0.8.12) and
$C_{k}(\cdot)$ is the $k$-th compound (0.8.1). Using the following properties:
$\operatorname{adj}_{0}(M)=\det(M)$ and $C_{0}(M)=1$, we have

$\det(M+\epsilon
B)=\epsilon^{0}\operatorname{tr}(\operatorname{adj}_{0}(M)C_{0}(B))+O(\epsilon)=\det(M)+O(\epsilon),$
(21)

completing the proof. ∎

## References

  * [1] Simon Barthelmé and Konstantin Usevich  “Spectral properties of kernel matrices in the flat limit”  In _SIAM Journal on Matrix Analysis and Applications_ 42.1, 2021, pp. 17–57  DOI: 10.1137/19M129677X
  * [2] D. Batenkov and Y. Yomdin  “On the Accuracy of Solving Confluent Prony Systems”  In _SIAM J. Appl. Math._ 73.1, 2013, pp. 134–154  DOI: 10.1137/110836584
  * [3] Dmitry Batenkov  “Decimated Generalized Prony Systems”, 2013  arXiv: http://arxiv.org/abs/1308.0753
  * [4] Dmitry Batenkov  “Stability and Super-Resolution of Generalized Spike Recovery”  In _Applied and Computational Harmonic Analysis_ 45.2, 2018, pp. 299–323  DOI: 10.1016/j.acha.2016.09.004
  * [5] Dmitry Batenkov, Laurent Demanet, Gil Goldman and Yosef Yomdin  “Conditioning of Partial Nonuniform Fourier Matrices with Clustered Nodes”  In _SIAM Journal on Matrix Analysis and Applications_ 44.1, 2020, pp. 199–220  DOI: 10/ggjwzb
  * [6] Dmitry Batenkov and Nuha Diab  “Super-Resolution of Generalized Spikes and Spectra of Confluent Vandermonde Matrices”  In _Applied and Computational Harmonic Analysis_ 65, 2023, pp. 181–208  DOI: 10.1016/j.acha.2023.03.002
  * [7] Dmitry Batenkov, Benedikt Diederichs, Gil Goldman and Yosef Yomdin  “The Spectral Properties of Vandermonde Matrices with Clustered Nodes”  In _Linear Algebra and its Applications_ 609, 2021, pp. 37–72  DOI: 10.1016/j.laa.2020.08.034
  * [8] Dmitry Batenkov, Gil Goldman and Yosef Yomdin  “Super-Resolution of near-Colliding Point Sources”  In _Information and Inference: A Journal of the IMA_ 10.2 Oxford Academic, 2021, pp. 515–572  DOI: 10.1093/imaiai/iaaa005
  * [9] M. Bertero and C. Mol  “III Super-Resolution by Data Inversion”  In _Progress in Optics_ 36, 1996, pp. 129–178
  * [10] K.. Chung and T.. Yao  “On Lattices Admitting Unique Lagrange Interpolations”  In _SIAM Journal on Numerical Analysis_ 14.4, 1977, pp. 735–743  DOI: 10.1137/0714050
  * [11] Annie Cuyt and Wen-shin Lee  “How to Get High Resolution Results from Sparse and Coarsely Sampled Data”  In _Applied and Computational Harmonic Analysis_ , 2018  DOI: 10/ggb5cv
  * [12] Laurent Demanet and Nam Nguyen  “The Recoverability Limit for Superresolution via Sparsity”  In _arXiv preprint arXiv:1502.01385_ , 2015  arXiv:1502.01385
  * [13] Benedikt Diederichs, Mihail N. Kolountzakis and Effie Papageorgiou  “How Many Fourier Coefficients Are Needed?”  In _Monatshefte für Mathematik_ , 2022  DOI: 10.1007/s00605-022-01792-0
  * [14] D.L. Donoho  “Superresolution via Sparsity Constraints”  In _SIAM Journal on Mathematical Analysis_ 23.5, 1992, pp. 1309–1331
  * [15] T.A. Driscoll and B. Fornberg  “Interpolation in the Limit of Increasingly Flat Radial Basis Functions”  In _Computers & Mathematics with Applications_ 43.3-5, 2002, pp. 413–422  DOI: 10.1016/S0898-1221(01)00295-4
  * [16] Hernán García, Camilo Hernández, Mauricio Junca and Mauricio Velasco  “Approximate Super-Resolution of Positive Measures in All Dimensions”  In _Applied and Computational Harmonic Analysis_ , 2020  DOI: 10/ggkz93
  * [17] Mariano Gasca and Thomas Sauer  “Polynomial Interpolation in Several Variables”  In _Advances in Computational Mathematics_ 12.4, 2000, pp. 377–410  DOI: 10.1023/A:1018981505752
  * [18] Gene H Golub and Charles F Van Loan  “Matrix Computations”  JHU Press, 2013
  * [19] Roger A Horn and Charles R Johnson  “Matrix analysis”  Cambridge university press, 2013
  * [20] Rami Katz, Nuha Diab and Dmitry Batenkov  “Decimated Prony’s Method for Stable Super-Resolution”  In _IEEE Signal Processing Letters_ 30, 2023, pp. 1467–1471  DOI: 10.1109/LSP.2023.3324553
  * [21] Rami Katz, Nuha Diab and Dmitry Batenkov  “On the Accuracy of Prony’s Method for Recovery of Exponential Sums with Closely Spaced Exponents”, 2024  DOI: 10.48550/arXiv.2302.05883
  * [22] Stefan Kunis and Dominik Nagel  “On the Smallest Singular Value of Multivariate Vandermonde Matrices with Clustered Nodes”  In _Linear Algebra and its Applications_ 604, 2020, pp. 1–20  DOI: 10.1016/j.laa.2020.06.003
  * [23] Stefan Kunis, Dominik Nagel and Anna Strotmann  “Multivariate Vandermonde Matrices with Separated Nodes on the Unit Circle Are Stable”  In _Applied and Computational Harmonic Analysis_ 58, 2022, pp. 50–59  DOI: 10.1016/j.acha.2022.01.001
  * [24] Weilin Li  “Nonharmonic Multivariate Fourier Transforms and Matrices: Condition Numbers and Hyperplane Geometry”, 2024
  * [25] Weilin Li and Wenjing Liao  “Stable super-resolution limit and smallest singular value of restricted Fourier matrices”  In _Applied and Computational Harmonic Analysis_ 51, 2020, pp. 118–156  DOI: 10.1016/j.acha.2020.10.004
  * [26] Weilin Li, Wenjing Liao and Albert Fannjiang  “Super-Resolution Limit of the ESPRIT Algorithm”  In _IEEE Transactions on Information Theory_ , 2020  DOI: 10/ggrnpw
  * [27] J. Lindberg  “Mathematical concepts of optical superresolution”  In _Journal of Optics_ 14, 2012
  * [28] Ping Liu and Hai Zhang  “A Mathematical Theory of Computational Resolution Limit in Multi-Dimensional Spaces”  In _Inverse Problems_ , 2021  DOI: 10.1088/1361-6420/ac245b
  * [29] Ping Liu and Hai Zhang  “A Mathematical Theory of the Computational Resolution Limit in One Dimension”  In _Applied and Computational Harmonic Analysis_ 56, 2022, pp. 402–446  DOI: 10.1016/j.acha.2021.09.002
  * [30] Charles A. Micchelli  “Interpolation of Scattered Data: Distance Matrices and Conditionally Positive Definite Functions”  In _Constructive Approximation_ 2.1, 1986, pp. 11–22  DOI: 10.1007/BF01893414
  * [31] R.. Nicolaides  “On a Class of Finite Elements Generated by Lagrange Interpolation”  In _SIAM Journal on Numerical Analysis_ 9.3, 1972, pp. 435–445  DOI: 10.1137/0709039
  * [32] Juan Manuel Peña and Tomas Sauer  “Stabilized Recovery and Model Reduction for Multivariate Exponential Polynomials”  In _Journal of Symbolic Computation_ , 2024, pp. 102313  DOI: 10.1016/j.jsc.2024.102313
  * [33] Markus Petz, Gerlind Plonka and Nadiia Derevianko  “Exact Reconstruction of Sparse Non-Harmonic Signals from Their Fourier Coefficients”  In _Sampling Theory, Signal Processing, and Data Analysis_ 19.1, 2021, pp. 7  DOI: 10.1007/s43670-021-00007-1
  * [34] Clarice Poon and Gabriel Peyré  “MultiDimensional Sparse Super-Resolution”  In _SIAM Journal on Mathematical Analysis_ 51.1, 2019, pp. 1–44  DOI: 10.1137/17M1147822
  * [35] Walter Rudin  “Principles of mathematical analysis”  McGraw-hill New York, 1964
  * [36] Souleymen Sahnoun, Konstantin Usevich and Pierre Comon  “Multidimensional ESPRIT for Damped and Undamped Signals: Algorithm, Computations, and Perturbation Analysis”  In _IEEE Transactions on Signal Processing_ 65.22, 2017, pp. 5897–5910  DOI: 10.1109/TSP.2017.2736512
  * [37] Tomas Sauer  “Polynomial Interpolation in Several Variables: Lattices, Differences, and Ideals”  In _Studies in Computational Mathematics_ 12 Elsevier, 2006, pp. 191–230  DOI: 10.1016/S1570-579X(06)80009-1
  * [38] Tomas Sauer  “Prony’s Method in Several Variables”  In _Numerische Mathematik_ 136.2, 2017, pp. 411–438  DOI: 10.1007/s00211-016-0844-8
  * [39] Tomas Sauer  “Prony’s Method in Several Variables: Symbolic Solutions by Universal Interpolation”  In _Journal of Symbolic Computation_ 84, 2018, pp. 95–112  DOI: 10.1016/j.jsc.2017.03.006
  * [40] Robert Schaback  “Multivariate interpolation by polynomials and radial basis functions”  In _Constructive approximation_ 21.3 New York, NY: Springer, 2005, pp. 293–317
  * [41] Andrew J Wathen and Shengxin Zhu  “On spectral distribution of kernel matrices related to radial basis functions”  In _Numerical Algorithms_ 70.4 Springer, 2015, pp. 709–726


# DCAlign v1.0: Aligning biological sequences using co-evolution models and
informed priors

Anna Paola Muntoni ${}^{\text{1,2}*}$, Andrea Pagnani ${}^{\text{1,2,3}}$
${}^{\text{\sf 1}}$Italian Institute for Genomic Medicine, IRCCS Candiolo,
SP-142, I-10060 Candiolo (TO), Italy
${}^{\text{\sf 2}}$Politecnico di Torino, Corso Duca degli Abruzzi, 24,
I-10129, Torino, Italy
${}^{\text{\sf 3}}$INFN, Sezione di Torino, Torino, Italy

###### Abstract

Summary: DCAlign is a new alignment method able to cope with the conservation
and the co-evolution signals that characterize the columns of multiple
sequence alignments of homologous sequences. However, the pre-processing steps
required to align a candidate sequence are computationally demanding. We show
in v1.0 how to dramatically reduce the overall computing time by including an
empirical prior over an informative set of variables mirroring the presence of
insertions and deletions.
Availability and implementation: DCAlign v1.0 is implemented in Julia and it
is fully available at https://github.com/infernet-h2020/DCAlign
Contact<EMAIL_ADDRESS>
Supplementary information: Supplementary data are available at Bioinformatics
online.

##  I Introduction

A common task in Bioinformatics is to cast evolutionary-related biological
sequences into a Multiple Sequence Alignment (MSA). The objective of this task
is to identify and align conserved regions of the sequences by maximizing the
similarity among the columns of the MSA. State-of-the-art alignment methods,
like HMMER for proteins [4], and Infernal [12] for RNAs, use hand-curated MSAs
of small representative subsets of sequences to be aligned (the so-called seed
alignments). Whereas for proteins, HMMER builds the Hidden Markov Model (HMM)
by using only the seed alignment, Infernal needs also secondary structure
information to generate a Covariance Model (CM). In both cases, HMM (for
proteins) or CM (for RNAs) are used to align query sequences. However,
homologous sequences show signals of correlated mutations (epistasis)
undetected by profile models.
Conservation and co-evolution signals are at the basis of Direct Coupling
Analysis (DCA)-based statistical models [9, 3]. Recently, these models have
been used to align biological sequences [10] and perform remote homology
search [17] by alignment of the sequences to a seed model, or by pairwise
alignments of seed models [15]. The method in [10], viz. DCAlign, returns the
ordered sub-sequence of a query unaligned sequence which maximizes an
objective function related to the DCA model of the seed. In this latter case,
standard DCA models fail to adequately describe the statistics of insertions
and gaps.
To alleviate this limitation, we added to the objective function gap and
insertion penalties learned from the seed alignment. While for the insertions,
the computational complexity is negligible, inferring gap penalties is a time-
consuming problem (see [10] and Supplementary text). Here, we treat penalties
in terms of informed priors computed from the seed sequences. The parameters
for gaps and insertions, extracted from the seed alignment, are determined in
an unsupervised manner. Finally, to further speed up the learning of the seed-
based objective function, we obtain the parameters of the DCA model using
pseudo-likelihood maximization [5] instead of Boltzmann Machine Learning [6,
11]. DCAlign v1.0, is a computational pipeline that allows for the computation
of the seed-model parameters in a few minutes, contrary to its original
implementation which required at least a day of computation in the best
scenario. The alignment problem is then solved approximately through a
message-passing algorithm (see Supplementary text).

##  II Methods

Figure 1: Schematic representation of the DCAlign v1.0 pipeline. From a
(given) hand-curated alignment (the seed, shown in the left panel), our
algorithm learns (i) a DCA model $\mathcal{H}$ exploiting the one-site and
two-site statistics of the seed (upper central box), and (ii) the gap and
insertion penalties by means of the empirical distribution of the pointer
differences $P\left(\Delta n_{ij}\right)$ for $i=1,\ldots,L$, and $j>i$
(bottom central box). The three sequences represent the three scenarios that
can occur between position $i$ and $j=i+1$: some insertion can appear, no
insertion and no gap is present, or $i+1$ contain a gap, so $\Delta
n_{i,i+1}=0$. For $j>i+1$ (lighter blue cases), both insertions and matched
symbols contribute to the computation of the $\Delta n_{i,j}$, while gaps do
not carry any contribution (see Fig. S1 for a more detailed example). The
alignment problem is then mapped into a constrained optimization problem over
the $\left(\boldsymbol{x},\boldsymbol{n}\right)$ variables. The constraints on
the variables and an example of alignment are shown in the right panel.

Our alignment algorithm estimates the optimal ordered sub-sequence compatible
with a DCA model and empirical knowledge of insertions and gaps of the seed.
Let $\boldsymbol{A}$ be an unaligned sequence of length $N$, and
$\boldsymbol{S}$ be its aligned counterpart of length $L$ (which is the length
of the seed MSA). We only consider the $L\leq N$ case. At each $i=1,\ldots,L$,
we define a Boolean variable $x_{i}\in\\{0,1\\}$ and a pointer
$n_{i}\in\\{0,\ldots,N+1\\}$. The variable $x_{i}$ indicates whether the
position $i$ is a gap ‘-‘ ($x_{i}=0$) or a match, i.e. a symbol in
$\boldsymbol{A}$. When $i$ is a match, $n_{i}$ identifies where $S_{i}$
matches $\boldsymbol{A}$, i.e. $S_{i}=A_{n_{i}}$; instead, for $x_{i}=0$, the
value of $n_{i}$ is used for keeping track of the last matched symbol in
$\boldsymbol{A}$. Let us define a pointer-difference variable as $\Delta
n_{i,j}=n_{j}-n_{i}$ for $i=1,\ldots,L$ and $j>i$. Each auxiliary variable
$\Delta n_{i,j}$ quantifies how many symbols of the unaligned sequence
$\boldsymbol{A}$ are present between two $i,j$ positions of the aligned
counterpart $\boldsymbol{S}$. If a configuration of the $\boldsymbol{n}$ is
given, the full set of the pointer differences reveal the presence of
insertions and gaps between any columns $i$ and $j$ of the alignment (see
Supplementary text).

###  II.1 Seed modeling

Together with a DCA model of the aligned seed (see Fig.1, central panel), for
every site $i$ (in red), we compute the $\Delta n_{i,j}$ for $j>i$ for all the
seed sequences, and we learn an empirical probability $P_{i,j}\left(\Delta
n_{i,j}\right)$ as shown in the bottom central panel of Fig.1 (this procedure
is computationally very fast). The color gradient is associated with the value
of $j$, the lighter the color, the larger is $j$. In Fig.1 (bottom central
panel) we consider as an example three sequences differing in the nature of
the $\Delta n_{i,i+1}$.

###  II.2 Alignment procedure

We can express the alignment problem in terms of the following optimization
problem:

$\boldsymbol{x},\boldsymbol{n}=\mathrm{argmax}_{\boldsymbol{\bar{x}},\boldsymbol{\bar{n}}}\frac{e^{-\beta\mathcal{H}\left(\boldsymbol{\bar{x}},\boldsymbol{\bar{n}}\right)}}{Z(\beta)}\prod_{i,j}P_{ij}^{\beta}\left(\boldsymbol{\bar{n}}\right),$
(1)

where $\mathcal{H}$ is the DCA model describing the seed (see Fig. 1, top
central panel), $Z$ is a normalization factor, and $\beta$ is a free parameter
whose relevance will be discussed below. The maximization only runs over the
feasible assignment of the variables, i.e. we impose that $n_{i+1}>n_{i}$ for
every column $i$. The informed prior will guide the optimization process
towards solutions that, among those that maximize the Boltzmann distribution
associated with $\mathcal{H}$, reproduce the statistics of the seed pointer
differences. Unfortunately, the problem thus stated is unfeasible as the
normalization function $Z$ cannot be efficiently computed. Similarly to the
first DCAlign version, we use an approximate message-passing algorithm coupled
with an annealing scheme over $\beta$ (i.e. we iteratively increase $\beta$)
to get the best alignment for the query sequence $\boldsymbol{A}$ (see
Supplementary text and Fig. S2).

##  III Results

We can classify the type of tests performed to assess the performance of our
computational strategy into three different categories:

  * •

Comparison with the previous implementation: As in [10], we compared our
results against HMMER, Infernal (the last algorithm only for RNA sequences) on
four Pfam (PF00035, PF00677, PF00684, PF00763), and Rfam (RF00059, RF00162,
RF00167, RF01734) families. A detailed description of the dataset is contained
in Tabs. S2-3. We utilized the following comparison metrics: (i) the positive
predictive value (PPV) of the DCA-based contact prediction [9, 3], (ii) the
proximity measures between the generated and the seed MSAs. As far as the
contact map prediction is concerned, we observe either a mild improvement or a
similar performance. With respect to the proximity measures, we notice a
negligible increase in the average distance between seed sequences and
generated alignments (see Figs. S3-6, and Tabs. S7-10).

  * •

Leave-one-out experiment: As a stress test for DCAlign v1.0 we also compared
our results to twenty-five ground-truth MSAs either extracted from benchmark
sets [2, 16, 7] or built from structural alignments [1] (see Tabs. S2, S4-5).
The numerical experiments consist of iteratively excluding one of the
sequences of the reference alignment and training HMM, CM, or DCAlign using
the remaining sequences. The excluded sequence is then aligned and
quantitatively compared to the ground truth (viz. the structural alignment, or
the benchmark sets). The emerging picture depends on the data type considered:
for benchmark sets all computational strategies seem to perform reasonably
well. In particular, HMMER (resp. Infernal) and our algorithm provide similar
outcomes for protein (resp. RNA) domains (see Figs. S7-10, Tabs. S11-12).
However, when we consider structural alignments as our reference ground truth,
our method significantly outperforms HMMER as shown in Figs. S11-12 and Tabs.
S13-14.

  * •

Divergent sequence alignment: Finally, to assess our algorithm’s remote
homology detection performance, we considered three RNA benchmark sets (the
seed of Rfam RF00162 [8], Twister type P1 [13], tRNA [14], see Tab. S6) from
[17]. Results suggest that Infernal is the best-performing method on two of
the three datasets, while our method achieves the best alignment for the tRNA
case. Note that Infernal is trained using secondary structure information that
our algorithm does not use. All results are presented in Fig. S13 and Tab.
S15.

From a computational efficiency point of view, the time needed to train the
algorithm is significantly smaller than both our old implementation and CM-
Infernal (see Supplementary text and Fig. S14). However, the time necessary to
align a sequence is equivalent compared to DCAlign, and probably to other
computational strategies taking into account epistasis [17, 15].

##  IV Conclusion

DCAlign v1.0 is a new implementation of the DCA-based alignment technique,
DCAlign, which conversely to the first implementation, allows for a fast
parametrization of the seed alignment. The new modeling significantly drops
the pre-processing time and guarantees a qualitatively equivalent alignment of
a set of target sequences.

## Acknowledgements

APM and AP acknowledge financial support from Marie Skłodowska-Curie, grant
agreement no. 734439(INFERNET). We also warmly thank Indaco Biazzo, Alfredo
Braunstein, Louise Budzynski and Luca Dall’Asta for interesting discussions.

## Data availability

DCAlign v1.0. is available at https://github.com/infernet-h2020/DCAlign

## References

  * Akdel et al. [2020] Akdel, M., Durairaj, J., de Ridder, D., and van Dijk, A. D. (2020).  Caretta - a multiple protein structure alignment and feature extraction suite.  Computational and Structural Biotechnology Journal, 18, 981–992.
  * Bahr et al. [2001] Bahr, A., Thompson, J. D., Thierry, J.-C., and Poch, O. (2001).  BAliBASE (Benchmark Alignment dataBASE): enhancements for repeats, transmembrane sequences and circular permutations.  Nucleic Acids Research, 29(1), 323–326.
  * Cocco et al. [2018] Cocco, S., Feinauer, C., Figliuzzi, M., Monasson, R., and Weigt, M. (2018).  Inverse statistical physics of protein sequences: a key issues review.  Reports on Progress in Physics, 81(3), 032601.
  * Eddy [2011] Eddy, S. R. (2011).  Accelerated Profile HMM Searches.  PLOS Computational Biology, 7(10), e1002195.  Publisher: Public Library of Science.
  * Ekeberg et al. [2013] Ekeberg, M., Lövkvist, C., Lan, Y., Weigt, M., and Aurell, E. (2013).  Improved contact prediction in proteins: Using pseudolikelihoods to infer Potts models.  Physical Review E, 87(1), 012707.  Publisher: American Physical Society.
  * Figliuzzi et al. [2018] Figliuzzi, M., Barrat-Charlaix, P., and Weigt, M. (2018).  How Pairwise Coevolutionary Models Capture the Collective Residue Variability in Proteins?  Molecular Biology and Evolution, 35(4), 1018–1027.
  * Freyhult et al. [2007] Freyhult, E. K., Bollback, J. P., and Gardner, P. P. (2007).  Exploring genomic dark matter: A critical assessment of the performance of homology search methods on noncoding RNA.  Genome Research, 17(1), 117–125.
  * Kalvari et al. [2020] Kalvari, I., Nawrocki, E. P., Ontiveros-Palacios, N., Argasinska, J., Lamkiewicz, K., Marz, M., Griffiths-Jones, S., Toffano-Nioche, C., Gautheret, D., Weinberg, Z., Rivas, E., Eddy, S. R., Finn, R., Bateman, A., and Petrov, A. I. (2020).  Rfam 14: expanded coverage of metagenomic, viral and microRNA families.  Nucleic Acids Research, 49(D1), D192–D200.
  * Morcos et al. [2011] Morcos, F., Pagnani, A., Lunt, B., Bertolino, A., Marks, D. S., Sander, C., Zecchina, R., Onuchic, J. N., Hwa, T., and Weigt, M. (2011).  Direct-coupling analysis of residue coevolution captures native contacts across many protein families.  Proceedings of the National Academy of Sciences, 108(49), E1293–E1301.  Publisher: Proceedings of the National Academy of Sciences.
  * Muntoni et al. [2020] Muntoni, A. P., Pagnani, A., Weigt, M., and Zamponi, F. (2020).  Aligning biological sequences by exploiting residue conservation and coevolution.  Physical Review E, 102(6), 062409.  Publisher: American Physical Society.
  * Muntoni et al. [2021] Muntoni, A. P., Pagnani, A., Weigt, M., and Zamponi, F. (2021).  adabmDCA: adaptive Boltzmann machine learning for biological sequences.  BMC Bioinformatics, 22(1), 528.
  * Nawrocki and Eddy [2013] Nawrocki, E. P. and Eddy, S. R. (2013).  Infernal 1.1: 100-fold faster RNA homology searches.  Bioinformatics, 29(22), 2933–2935.
  * Roth et al. [2014] Roth, A., Weinberg, Z., Chen, A. G. Y., Kim, P. B., Ames, T. D., and Breaker, R. R. (2014).  A widespread self-cleaving ribozyme class is revealed by bioinformatics.  Nature Chemical Biology, 10(1), 56–60.  Number: 1 Publisher: Nature Publishing Group.
  * Sprinzl et al. [1998] Sprinzl, M., Horn, C., Brown, M., Ioudovitch, A., and Steinberg, S. (1998).  Compilation of tRNA sequences and sequences of tRNA genes.  Nucleic Acids Research, 26(1), 148–153.
  * Talibart and Coste [2021] Talibart, H. and Coste, F. (2021).  PPalign: optimal alignment of Potts models representing proteins with direct coupling information.  BMC Bioinformatics, 22(1).
  * Thompson et al. [2005] Thompson, J. D., Koehl, P., Ripp, R., and Poch, O. (2005).  BAliBASE 3.0: Latest developments of the multiple sequence alignment benchmark.  Proteins: Structure, Function, and Bioinformatics, 61(1), 127–136.
  * Wilburn and Eddy [2020] Wilburn, G. W. and Eddy, S. R. (2020).  Remote homology search with hidden Potts models.  PLOS Computational Biology, 16(11), e1008085.  Publisher: Public Library of Science.