# Fast and direct inversion methods for the multivariate nonequispaced fast
Fourier transform

Melanie Kircheis111Corresponding author<EMAIL_ADDRESS>Chemnitz University of Technology, Faculty of Mathematics, D–09107 Chemnitz,
Germany Daniel<EMAIL_ADDRESS>Chemnitz University of
Technology, Faculty of Mathematics, D–09107 Chemnitz, Germany

###### Abstract

The well-known discrete Fourier transform (DFT) can easily be generalized to
arbitrary points in the spatial domain. The fast procedure for this
generalization is referred to as nonequispaced fast Fourier transform (NFFT).
Various applications such as MRI, solution of PDEs, etc., are interested in
the inverse problem, i. e., computing Fourier coefficients from given
nonequispaced data. In this paper we survey different kinds of approaches to
tackle this problem. In contrast to iterative procedures, where multiple
iteration steps are needed for computing a solution, we focus especially on
so-called direct inversion methods. We review density compensation techniques
and introduce a new scheme that leads to an exact reconstruction for
trigonometric polynomials. In addition, we consider a matrix optimization
approach using Frobenius norm minimization to obtain an inverse NFFT.

_Key words_ : inverse nonequispaced fast Fourier transform, nonuniform fast
Fourier transform, direct inversion, density compensation, matrix
optimization, iNFFT, NFFT, NUFFT

AMS _Subject Classifications_ :  65Txx, 65T50, 65F05.

##  1 Introduction

The NFFT, short hand for nonequispaced fast Fourier transform or nonuniform
fast Fourier transform (NUFFT), respectively, is a fast algorithm to evaluate
a trigonometric polynomial

$f(\boldsymbol{x})=\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\hat{f}_{\boldsymbol{k}}\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}}$
(1.1)

with given Fourier coefficients $\hat{f}_{\boldsymbol{k}}\in\mathbb{C}$,
$\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}$, at nonequispaced points
$\boldsymbol{x}_{j}\in\left[-\frac{1}{2},\frac{1}{2}\right)^{d}$,
$j=1,\dots,N$, $N\in\mathbb{N}$, where
$\mathcal{I}_{\boldsymbol{M}}\coloneqq\mathbb{Z}^{d}\cap\left[-\tfrac{M}{2},\tfrac{M}{2}\right)^{d}$
with $|\mathcal{I}_{\boldsymbol{M}}|=M^{d}$. In case we are given equispaced
points $\boldsymbol{x}_{j}$ and $|\mathcal{I}_{\boldsymbol{M}}|=N$, this
evaluation can be realized by means of the well-known fast Fourier transform
(FFT); an algorithm that is invertible. However, various applications such as
magnetic resonance imaging (MRI), cf. [15, 19], solution of PDEs, cf. [21],
etc., need to perform an inverse nonequispaced fast Fourier transform (iNFFT),
i. e., compute the Fourier coefficients $\hat{f}_{\boldsymbol{k}}$ from given
function evaluations $f(\boldsymbol{x}_{j})$ of the trigonometric polynomial
(1.1). Hence, we are interested in an inversion also for nonequispaced data.

In general, the number $N$ of points $\boldsymbol{x}_{j}$ is independent from
the number $|\mathcal{I}_{\boldsymbol{M}}|$ of Fourier coefficients
$\hat{f}_{\boldsymbol{k}}$ and hence the nonequispaced Fourier matrix

$\boldsymbol{A}\coloneqq\left(\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}}\right)_{j=1,\,\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}^{N}\
\in\mathbb{C}^{N\times|\mathcal{I}_{\boldsymbol{M}}|},$

which we would have to invert, is rectangular in most cases. Considering the
corresponding linear system
$\boldsymbol{A}\boldsymbol{\hat{f}}=\boldsymbol{f}$ with
$\boldsymbol{f}\coloneqq\left(f(\boldsymbol{x}_{j})\right)_{j=1}^{N}$ and
$\boldsymbol{\hat{f}}\coloneqq(\hat{f}_{\boldsymbol{k}})_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}$,
this can either be overdetermined, if $|\mathcal{I}_{\boldsymbol{M}}|\leq N$,
or underdetermined, if $|\mathcal{I}_{\boldsymbol{M}}|>N$. Generally, this
forces us to look for a pseudoinverse solution. Moreover, we also require that
the nonequispaced Fourier matrix $\boldsymbol{A}$ has full rank. Eigenvalue
estimates in [22, 6, 9, 45, 43] indeed confirm that this condition is
satisfied for sufficiently nice sampling sets.

In literature a variety of approaches for an inverse NFFT (iNFFT) can be
found. This is why we give a short overview.

##### Iterative inversion methods

We start surveying the iterative methods. For the one-dimensional setting
$d=1$ with $|\mathcal{I}_{\boldsymbol{M}}|=N$ an algorithm was published in
[60], which is specially designed for jittered equispaced points and is based
on the conjugate gradient (CG) method in connection with low rank
approximation. An approach for the overdetermined case
$|\mathcal{I}_{\boldsymbol{M}}|\leq N$ can be found in [22], where the
Toeplitz structure of the matrix product
$\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A}$ with a diagonal matrix
$\boldsymbol{W}\coloneqq\operatorname*{diag}(w_{j})_{j=1}^{N}$ of Voronoi
weights is used to compute the solution iteratively by means of the CG
algorithm.

For higher dimensional problems with $d\geq 1$ there are several approaches
that compute a least squares approximation to the linear system
$\boldsymbol{A}\boldsymbol{\hat{f}}=\boldsymbol{f}$. In the overdetermined
case $|\mathcal{I}_{\boldsymbol{M}}|\leq N$, the given data can typically only
be approximated up to a residual
$\boldsymbol{r}\coloneqq\boldsymbol{A}\boldsymbol{\hat{f}}-\boldsymbol{f}$.
Therefore, the weighted least squares problem

$\displaystyle\underset{\boldsymbol{\hat{f}}\in\mathbb{C}^{|\mathcal{I}_{\boldsymbol{M}}|}}{\text{Minimize
}}\,\,\sum_{j=1}^{N}w_{j}\,\Bigg{|}\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\hat{f}_{\boldsymbol{k}}\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}}–f(\boldsymbol{x}_{j})\Bigg{|}^{2}$

is considered, which is equivalent to solving the weighted normal equations of
first kind
$\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A}\boldsymbol{\hat{f}}=\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{f}$
with the diagonal matrix
$\boldsymbol{W}\coloneqq\operatorname*{diag}(w_{j})_{j=1}^{N}$ of weights in
time domain. In [69, 24, 41] these normal equations are solved iteratively by
means of CG using the NFFT to realize fast matrix-vector multiplications
involving $\boldsymbol{A}$, whereas in [56] a fast convolution is used.

In the consistent underdetermined case $|\mathcal{I}_{\boldsymbol{M}}|>N$ the
data can be interpolated exactly and therefore one can choose a specific
solution, e. g. the one that solves the constrained minimization problem

$\displaystyle\underset{\boldsymbol{\hat{f}}\in\mathbb{C}^{|\mathcal{I}_{\boldsymbol{M}}|}}{\text{Minimize
}}\,\,\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\frac{|\hat{f}_{\boldsymbol{k}}|^{2}}{\hat{w}_{\boldsymbol{k}}}\quad\text{
subject to }\quad\boldsymbol{A}\boldsymbol{\hat{f}}=\boldsymbol{f}.$

This interpolation problem is equivalent to the weighted normal equations of
second kind
$\boldsymbol{A}\boldsymbol{\hat{W}}\boldsymbol{A}^{*}\boldsymbol{y}=\boldsymbol{f}$,
$\boldsymbol{\hat{f}}=\boldsymbol{\hat{W}}\boldsymbol{A}^{*}\boldsymbol{y}$
with the diagonal matrix
$\boldsymbol{\hat{W}}\coloneqq\operatorname*{diag}(\hat{w}_{\boldsymbol{k}})_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}$
of weights in frequency domain. In [44] the CG method was used in connection
with the NFFT to iteratively compute a solution to this problem, see also [53,
Section 7.6.2].

##### Regularization methods

Moreover, there also exist several regularization techniques for the
multidimensional setting $d\geq 1$. For example, [65, 71, 1] all solve the
$\ell_{1}$-regularized problem

$\displaystyle\underset{\boldsymbol{\hat{f}}\in\mathbb{C}^{|\mathcal{I}_{\boldsymbol{M}}|}}{\text{Minimize
}}\,\tfrac{1}{2}\|\boldsymbol{A}\boldsymbol{\hat{f}}-\boldsymbol{f}\|_{2}^{2}+\lambda\|\mathcal{L}^{m}\boldsymbol{\hat{f}}\|_{1}$

with regularization parameter $\lambda>0$ and the $m$-th order polynomial
annihilation operator
$\mathcal{L}^{m}\in\mathbb{R}^{N\times|\mathcal{I}_{\boldsymbol{M}}|}$ as
sparsifying transform, see [2]. Based on this, weighted $\ell_{p}$-schemes

$\displaystyle\underset{\boldsymbol{\hat{f}}\in\mathbb{C}^{|\mathcal{I}_{\boldsymbol{M}}|}}{\text{Minimize
}}\,\tfrac{1}{2}\|\boldsymbol{A}\boldsymbol{\hat{f}}-\boldsymbol{f}\|_{2}^{2}+\tfrac{1}{p}\|\boldsymbol{W}\mathcal{L}^{m}\boldsymbol{\hat{f}}\|_{p}^{p}$

were introduced in [10, 11, 14, 48], which are designed to reduce the penalty
at locations where $\mathcal{L}^{m}\boldsymbol{\hat{f}}$ is nonzero. For
instance, [13, 61] each state a two step method, that firstly uses edge
detection to create a mask, i. e., a weighting matrix which indicates where
non-zero entries are expected in the TV domain, and then targets weighted
$\ell_{2}$-norm TV regularization appropriately to smooth regions of the
function in a second minimization step.

##### Direct inversion methods

In contrast to these iterated procedures, there are also so-called direct
methods that do not require multiple iteration steps. Already in [18] a direct
method was explained for the setting $d=1$ and
$|\mathcal{I}_{\boldsymbol{M}}|=N$ which uses Lagrange interpolation in
combination with fast multiple methods. Based on this, further methods were
deduced for the same setting, which also use Lagrange interpolation, but
additionally incorporate an imaginary shift in [63], or utilize fast summation
in [37] for the fast evaluation of occurring sums, see also [53, Section
7.6.1].

In the overdetermined setting $|\mathcal{I}_{\boldsymbol{M}}|\leq N$ another
approach for computing an inverse NFFT can be obtained by using the fact that
$\boldsymbol{A}^{*}\boldsymbol{A}$ is of Toeplitz structure. To this end, the
Gohberg-Semencul formula, see [32], can be used to solve the normal equations
$\boldsymbol{A}^{*}\boldsymbol{A}\boldsymbol{\hat{f}}=\boldsymbol{A}^{*}\boldsymbol{f}$
exactly by analogy with [5]. Here the computation of the components of the
Gohberg-Semencul formula can be viewed as a precomputational step. In
addition, also a frame-theoretical approach is known from [26], which provides
a link between the adjoint NFFT and frame approximation, and could therefore
be seen as a way to invert the NFFT. Note that the method in [26] is based
only on optimizing a diagonal matrix (the matrix $\boldsymbol{D}$ defined in
(2.13)), whereas in [37] similar ideas were used to modify a sparse matrix
(the matrix $\boldsymbol{B}$ defined in (2.15)).

For the multidimensional setting $d>1$ several methods have been developed
that are tailored to the special structure of the linogram or pseudo-polar
grid, respectively, see Figure 2(b), such that the inversion involves only
one-dimensional FFTs and interpolations. On the one hand, in [4] a least
squares solution is computed iteratively by using a fast multiplication
technique with the matrix $\boldsymbol{A}$, which can be derived in case of
the linogram grid. On the other hand, [3] is based on a fast resampling
strategy, where a first step resamples the linogram grid onto a Cartesian
grid, and the second phase recovers the image from these Cartesian samples.
However, these techniques are exclusively applicable for the special case of
the linogram grid, see Figure 2(b), or the polar grid by another resampling,
cf. Figure 2(a). Since we are interested in more generally applicable methods,
a brief introduction to direct inversion for general sampling patterns can be
found in [38].

##### Current approach

In this paper we focus on direct inversion methods that are applicable for
general sampling patterns and present new methods for the computation of an
iNFFT. Note that direct method in the context of the linear system
$\boldsymbol{A}\boldsymbol{\hat{f}}=\boldsymbol{f}$ means, that for a fixed
set of points $\boldsymbol{x}_{j}$, $j=1,\dots,N$, the reconstruction of
$\boldsymbol{\hat{f}}$ from given $\boldsymbol{f}$ can be realized with the
same number of arithmetic operations as a single application of an adjoint
NFFT (see Algorithm 2.4). To achieve this, a certain precomputational step is
compulsory, since the adjoint NFFT does not yield an inversion of the NFFT per
se, cf. (3.3). Although this precomputations might be rather costly, they need
to be done only once for a given set of points $\boldsymbol{x}_{j}$,
$j=1,\dots,N$. Therefore, direct methods are especially beneficial in case of
fixed points for several measurement vectors $\boldsymbol{f}$.

For this reason, the current paper is concerned with two different approaches
of this kind. Firstly, we consider the very well known approach of so-called
sampling density compensation, which can be written as
$\boldsymbol{\hat{f}}\approx\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{f}$
with a diagonal matrix
$\boldsymbol{W}\coloneqq\operatorname*{diag}(w_{j})_{j=1}^{N}$ of weights. The
already mentioned precomputations then consist of computing suitable weights
$w_{j}$, while the actual reconstruction step includes only one adjoint NFFT
applied to the scaled coefficient vector $\boldsymbol{W}\boldsymbol{f}$. In
this paper we examine several existing approaches for computing the weights
$w_{j}$ and introduce a new method, such that the iNFFT is exact for all
trigonometric polynomials (1.1) of degree $\boldsymbol{M}$.

As a second part, we reconsider and enhance our approach introduced in [38].
Here the idea is using the matrix representation
$\boldsymbol{A}\approx\boldsymbol{B}\boldsymbol{F}\boldsymbol{D}$ of the NFFT
to modify one of the matrix factors, such that an inversion is given by
$\boldsymbol{\hat{f}}\approx\boldsymbol{D}^{*}\boldsymbol{F}^{*}\boldsymbol{B}_{\mathrm{opt}}^{*}\boldsymbol{f}$.
Then the precomputational step consists of optimizing the matrix
$\boldsymbol{B}$, while the actual reconstruction step includes only one
modified adjoint NFFT applied to the coefficient vector $\boldsymbol{f}$.

##### Outline of this paper

This paper is organized as follows. In Section 2 we introduce the already
mentioned algorithm, the NFFT, as well as its adjoint version, the adjoint
NFFT. Secondly, in Section 3 we introduce the inversion problem and deal with
direct methods using so-called sampling density compensation. We start our
investigations with trigonometric polynomials in Section 3.1. Here the main
formula (3.16) yields exact reconstruction for all trigonometric polynomials
of degree $\boldsymbol{M}$. We also discuss practical computation schemes for
the overdetermined as well as the underdetermined setting. Subsequently, in
Section 3.2 we go on to bandlimited functions and show that the same numerical
procedures as in Section 3.1 can be used in this setting as well. Section 3.3
then summarizes the previous findings by presenting a general error bound on
density compensation factors computed by means of (3.16) in Theorem 3.14. In
addition, this also yields an estimate on the condition number of a specific
matrix product, as shown in Theorem 3.15. In Section 3.4 we have a look at
certain approaches from literature and their connection among each other as
well as to the method presented in Section 3.1. Afterwards, we examine another
direct inversion method in Section 4, where we aim to modify the adjoint NFFT
by means of matrix optimization such that this yields an iNFFT. Finally, in
Section 5 we show some numerical examples to investigate the accuracy of our
approaches.

##  2 Nonequispaced fast Fourier transform

Let

$\mathbb{T}^{d}\coloneqq\mathbb{R}^{d}\setminus\mathbb{Z}^{d}\cong\left[-\tfrac{1}{2},\tfrac{1}{2}\right)^{d}=\left\\{\boldsymbol{x}\in\mathbb{R}^{d}\colon-\tfrac{1}{2}\leq
x_{t}<\tfrac{1}{2},\,t=1,\dots,d\right\\}$

denote the $d$-dimensional torus with $d\in\mathbb{N}$. For
$\boldsymbol{M}\coloneqq(M,\dots,M)^{T}$, $M\in 2\mathbb{N}$, we define the
multi-index set

$\mathcal{I}_{\boldsymbol{M}}\coloneqq\mathbb{Z}^{d}\cap\left[-\tfrac{M}{2},\tfrac{M}{2}\right)^{d}=\left\\{\boldsymbol{k}\in\mathbb{Z}^{d}\colon-\tfrac{M}{2}\leq
k_{t}<\tfrac{M}{2},\,t=1,\dots,d\right\\}$

with cardinality $|\mathcal{I}_{\boldsymbol{M}}|=M^{d}$. The inner product of
two vectors shall be defined as usual as
$\boldsymbol{k}\boldsymbol{x}\coloneqq k_{1}x_{1}+\dots+k_{d}x_{d}$.
Additionally, we define the componentwise product as
$\boldsymbol{x}\odot\boldsymbol{y}\coloneqq\left(x_{1}y_{1},\dots,x_{d}y_{d}\right)^{T}$,
the all ones vector $\boldsymbol{1}_{d}\coloneqq(1,\dots,1)^{T}$ and the
reciprocal of a vector $\boldsymbol{x}$ with nonzero components shall be given
by $\boldsymbol{x}^{-1}\coloneqq\left(x_{1}^{-1},\dots,x_{d}^{-1}\right)^{T}$.

We consider the Hilbert space $L_{2}(\mathbb{T}^{d})$ of all $1$-periodic,
complex-valued functions, which possesses the orthonormal basis
$\\{\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}}\colon\boldsymbol{k}\in\mathbb{Z}^{d}\\}$.
It it known that every function $f\in L_{2}(\mathbb{T}^{d})$ is uniquely
representable in the form

$\displaystyle
f(\boldsymbol{x})=\sum_{\boldsymbol{k}\in\mathbb{Z}^{d}}c_{\boldsymbol{k}}(f)\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}}$
(2.1)

with the coefficients

$\displaystyle
c_{\boldsymbol{k}}(f)\coloneqq\int_{\mathbb{T}^{d}}f(\boldsymbol{x})\,\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}}\,\mathrm{d}\boldsymbol{x},\quad\boldsymbol{k}\in\mathbb{Z}^{d},$
(2.2)

where the sum in (2.1) converges to $f$ in the $L_{2}(\mathbb{T}^{d})$-norm,
cf. [53, Thm. 4.5]. A series of the form (2.1) is called _Fourier series_ with
the _Fourier coefficients_ (2.2). Numerically, the Fourier coefficients (2.2)
are approximated on the uniform grid
$\\{\boldsymbol{M}^{-1}\odot\boldsymbol{\ell},\,\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}}\\}$
by the trapezoidal rule for numerical integration as

$\displaystyle
c_{\boldsymbol{k}}(f)\approx\frac{1}{|\mathcal{I}_{\boldsymbol{M}}|}\sum_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}}}f(\boldsymbol{M}^{-1}\odot\boldsymbol{\ell})\,\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{k}(\boldsymbol{M}^{-1}\odot\boldsymbol{\ell})},\quad\boldsymbol{k}\in\mathbb{Z}^{d},$
(2.3)

which is an acceptable approximation for
$\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}$, see e. g. [53, p. 214]. The
fast evaluation of (2.3) can then be realized by means of the _fast Fourier
transform (FFT)_. Moreover, it is know that this transformation is invertible
and that the inverse problem of computing

$f(\boldsymbol{M}^{-1}\odot\boldsymbol{\ell})=\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\hat{f}_{\boldsymbol{k}}\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}(\boldsymbol{M}^{-1}\odot\boldsymbol{\ell})},\quad\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}},$
(2.4)

with $\hat{f}_{\boldsymbol{k}}\approx c_{\boldsymbol{k}}(f)$,
$\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}$, can be realized by means of
an _inverse fast Fourier transform (iFFT)_ , which is basically the same
algorithm except for some reordering and scaling, cf. [53, Lem. 3.17].

Now suppose we are given nonequispaced points
$\boldsymbol{x}_{j}\in\mathbb{T}^{d}$, $j=1,\dots,N$, instead. Then, we
consider the computation of the sums

$f_{j}\coloneqq
f(\boldsymbol{x}_{j})=\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\hat{f}_{\boldsymbol{k}}\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}},\quad
j=1,\dots,N,$ (2.5)

for given $\hat{f}_{\boldsymbol{k}}\in\mathbb{C}$,
$\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}$, as well as the adjoint
problem of computing

$h_{\boldsymbol{k}}=\sum_{j=1}^{N}f_{j}\,\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}},\quad\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}},$
(2.6)

for given $f_{j}\in\mathbb{C}$, $j=1,\dots,N$. By defining the nonequispaced
Fourier matrix

$\boldsymbol{A}=\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{M}}|}\coloneqq\left(\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}}\right)_{j=1,\,\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}^{N}\
\in\mathbb{C}^{N\times|\mathcal{I}_{\boldsymbol{M}}|},$ (2.7)

as well as the vectors $\boldsymbol{f}\coloneqq\left(f_{j}\right)_{j=1}^{N}$,
$\boldsymbol{\hat{f}}\coloneqq(\hat{f}_{\boldsymbol{k}})_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}$
and
$\boldsymbol{h}\coloneqq(h_{\boldsymbol{k}})_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}$,
the computation of sums of the form (2.5) and (2.6) can be written as
$\boldsymbol{f}=\boldsymbol{A}\boldsymbol{\hat{f}}$ and
$\boldsymbol{h}=\boldsymbol{A}^{*}\boldsymbol{f}$, where
$\boldsymbol{A}^{*}\coloneqq\overline{\boldsymbol{A}}^{T}$ denotes the adjoint
matrix of $\boldsymbol{A}$.

Since the naive computation of (2.5) and (2.6) is of complexity
$\mathcal{O}(N\cdot|\mathcal{I}_{\boldsymbol{M}}|)$, a fast approximate
algorithm, the so-called _nonequispaced fast Fourier transform (NFFT)_ , is
briefly explained below. For more information see [17, 7, 66, 29, 36] or [53,
pp. 377-381].

###  2.1 The NFFT

We firstly restrict our attention to problem (2.5), which is equivalent to the
evaluation of a trigonometric polynomial

$f(\boldsymbol{x})=\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\hat{f}_{\boldsymbol{k}}\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}}$
(2.8)

with given $\hat{f}_{\boldsymbol{k}}\in\mathbb{C}$,
$\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}$, at given nonequispaced points
$\boldsymbol{x}_{j}\in\mathbb{T}^{d}$, $j=1,\dots,N$. Let $w\in
L_{2}(\mathbb{R}^{d})\cap L_{1}(\mathbb{R}^{d})$ be a so-called _window
function_ , which is well localized in space and frequency domain. Now we
define the 1-periodic function
$\tilde{w}(\boldsymbol{x})\coloneqq\sum_{\boldsymbol{r}\in\mathbb{Z}^{d}}w(\boldsymbol{x}+\boldsymbol{r})$
with absolute convergent Fourier series. As a consequence, the Fourier
coefficients of the periodization $\tilde{w}$ have the form

$c_{\boldsymbol{k}}(\tilde{w})=\int_{\mathbb{T}^{d}}\tilde{w}(\boldsymbol{x})\,\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}}\,\mathrm{d}\boldsymbol{x}=\int_{\mathbb{R}^{d}}w(\boldsymbol{x})\,\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}}\,\mathrm{d}\boldsymbol{x}\eqqcolon\hat{w}(\boldsymbol{k}).$

For a given oversampling factor $\sigma\geq 1$, we define $2\mathbb{N}\ni
M_{\sigma}\coloneqq 2\lceil\lceil\sigma M\rceil/2\rceil$ as well as
$\boldsymbol{M}_{\boldsymbol{\sigma}}\coloneqq
M_{\sigma}\cdot\boldsymbol{1}_{d}$, and approximate $f$ by a linear
combination of translates of the periodized window function, i. e.,

$f(\boldsymbol{x})\approx
s_{1}(\boldsymbol{x})\coloneqq\sum_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}}g_{\boldsymbol{\ell}}\,\tilde{w}\left(\boldsymbol{x}-\boldsymbol{M}_{\boldsymbol{\sigma}}^{-1}\odot\boldsymbol{\ell}\right),$
(2.9)

where $g_{\boldsymbol{\ell}}\in\mathbb{C}$,
$\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}$, are
coefficients to be determined such that (2.9) yields a good approximation. By
means of the convolution theorem (see [53, Lem. 4.1]), the approximant
$s_{1}\in L_{2}(\mathbb{T}^{d})$ in (2.9) can be represented as

$\displaystyle s_{1}(\boldsymbol{x})$
$\displaystyle=\sum_{\boldsymbol{k}\in\mathbb{Z}^{d}}c_{\boldsymbol{k}}(s_{1})\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}}$
$\displaystyle=\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\hat{g}_{\boldsymbol{k}}\;c_{\boldsymbol{k}}(\tilde{w})\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}}+\sum_{\boldsymbol{r}\in\mathbb{Z}^{d}\setminus\\{\boldsymbol{0}\\}}\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\hat{g}_{\boldsymbol{k}}\;c_{\boldsymbol{k}+\boldsymbol{M}_{\boldsymbol{\sigma}}\odot\,\boldsymbol{r}}(\tilde{w})\,\mathrm{e}^{2\pi\mathrm{i}(\boldsymbol{k}+\boldsymbol{M}_{\boldsymbol{\sigma}}\odot\,\boldsymbol{r})\boldsymbol{x}},$
(2.10)

where the discrete Fourier transform of the coefficients
$g_{\boldsymbol{\ell}}$ is defined by

$\hat{g}_{\boldsymbol{k}}\coloneqq\sum_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}}g_{\boldsymbol{\ell}}\,\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{k}(\boldsymbol{M}_{\boldsymbol{\sigma}}^{-1}\odot\,\boldsymbol{\ell})},\quad\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}.$
(2.11)

Comparing (2.5) and (2.1) then yields

$\hat{g}_{\boldsymbol{k}}=\left\\{\begin{array}[]{cl}\dfrac{\hat{f}_{\boldsymbol{k}}}{\hat{w}(\boldsymbol{k})}&:\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}},\\\
0&:\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}\setminus\mathcal{I}_{\boldsymbol{M}}.\end{array}\right.$

Consequently, the coefficients $g_{\boldsymbol{\ell}}$ in (2.9) can be
obtained by inverting (2.11), i. e., by the application of an iFFT.

Furthermore, we assume that $w$ is well localized such that it is small
outside the square
$\left[-\nicefrac{{m}}{{M_{\sigma}}},\nicefrac{{m}}{{M_{\sigma}}}\right]^{d}$
with truncation parameter $m\ll M_{\sigma}$. In this case, $w$ can be
approximated by the compactly supported function

$w_{m}(\boldsymbol{x})\coloneqq\left\\{\begin{array}[]{cl}w(\boldsymbol{x})&:\boldsymbol{x}\in\left[-\frac{m}{M_{\sigma}},\frac{m}{M_{\sigma}}\right]^{d},\\\
0&:\text{otherwise}.\end{array}\right.$

Thereby, we approximate $s_{1}$ by the short sums

$\displaystyle f(\boldsymbol{x}_{j})\approx s_{1}(\boldsymbol{x}_{j})\approx
s(\boldsymbol{x}_{j})$
$\displaystyle\coloneqq\sum_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}}g_{\boldsymbol{\ell}}\,\tilde{w}_{m}\left(\boldsymbol{x}_{j}-\boldsymbol{M}_{\boldsymbol{\sigma}}^{-1}\odot\boldsymbol{\ell}\right)$
$\displaystyle=\sum_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}},m}(\boldsymbol{x}_{j})}g_{\boldsymbol{\ell}}\,\tilde{w}_{m}\left(\boldsymbol{x}_{j}-\boldsymbol{M}_{\boldsymbol{\sigma}}^{-1}\odot\boldsymbol{\ell}\right),$

where the index set

$\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}},m}(\boldsymbol{x}_{j})\coloneqq\left\\{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}\colon\exists\boldsymbol{z}\in\mathbb{Z}^{d}\text{
with
}-m\cdot\boldsymbol{1}_{d}\leq\boldsymbol{M}_{\boldsymbol{\sigma}}\odot\boldsymbol{x}_{j}-\boldsymbol{\ell}+\boldsymbol{z}\leq
m\cdot\boldsymbol{1}_{d}\right\\}$ (2.12)

contains at most $(2m+1)^{d}$ entries for each fixed $\boldsymbol{x}_{j}$.
Thus, the obtained algorithm can be summarized as follows.

######  Algorithm 2.1 (NFFT).

For $d,N\in\mathbb{N}$ let
$\boldsymbol{x}_{j}\in\mathbb{T}^{d},\,j=1,\dots,N,$ be given points as well
as  $\hat{f}_{\boldsymbol{k}}\in\mathbb{C}$,
$\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}$, given Fourier coefficients.
Furthermore, we are given the oversampling factor $\sigma\geq 1$,
$2\mathbb{N}\ni M_{\sigma}\coloneqq 2\lceil\lceil\sigma M\rceil/2\rceil$,
$\boldsymbol{M}_{\boldsymbol{\sigma}}\coloneqq
M_{\sigma}\cdot\boldsymbol{1}_{d}$, as well as the window function $w$, the
truncated function $w_{m}$ with truncation parameter $m\ll M_{\sigma},$ and
their 1-periodic versions $\tilde{w}$ and $\tilde{w}_{m}$.

  1. 1.

Set

$\hat{g}_{\boldsymbol{k}}\coloneqq\left\\{\begin{array}[]{cl}\frac{\hat{f}_{\boldsymbol{k}}}{\hat{w}(\boldsymbol{k})}&:\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}},\\\
0&:\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}\setminus\mathcal{I}_{\boldsymbol{M}}.\end{array}\right.$

$\mathcal{O}(|\mathcal{I}_{\boldsymbol{M}}|)$

  2. 2.

Compute

$g_{\boldsymbol{\ell}}\coloneqq\frac{1}{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|}\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\hat{g}_{\boldsymbol{k}}\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}(\boldsymbol{M}_{\boldsymbol{\sigma}}^{-1}\odot\,\boldsymbol{\ell})},\quad\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}},$

by means of a $d$-variate iFFT.
$\mathcal{O}(|\mathcal{I}_{\boldsymbol{M}}|\log(|\mathcal{I}_{\boldsymbol{M}}|))$

  3. 3.

Compute the short sums

$\tilde{f}_{j}\coloneqq\sum_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}},m}(\boldsymbol{x}_{j})}g_{\boldsymbol{\ell}}\,\tilde{w}_{m}\left(\boldsymbol{x}_{j}-\boldsymbol{M}_{\boldsymbol{\sigma}}^{-1}\odot\boldsymbol{\ell}\right),\quad
j=1,\dots,N.$

$\mathcal{O}(N)$

Output: $\tilde{f}_{j}\approx f_{j}$, $j=1,\dots,N$, cf. (2.5). Complexity:
$\mathcal{O}(|\mathcal{I}_{\boldsymbol{M}}|\log(|\mathcal{I}_{\boldsymbol{M}}|)+N)$

######  Remark 2.2.

Suitable window functions can be found e. g. in [17, 7, 16, 25, 29, 36, 54].

Next we give the matrix-vector representation of the NFFT. To this end, we
define

  * •

the diagonal matrix

$\boldsymbol{D}\coloneqq\text{diag}\left(\frac{1}{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|\cdot\hat{w}(\boldsymbol{k})}\right)_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\
\in\mathbb{C}^{|\mathcal{I}_{\boldsymbol{M}}|\times|\mathcal{I}_{\boldsymbol{M}}|},$
(2.13)

  * •

the truncated Fourier matrix

$\boldsymbol{F}\coloneqq\left(\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}(\boldsymbol{M}_{\boldsymbol{\sigma}}^{-1}\odot\,\boldsymbol{\ell})}\right)_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}},\,\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\
\in\mathbb{C}^{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|\times|\mathcal{I}_{\boldsymbol{M}}|},$
(2.14)

  * •

and the sparse matrix

$\boldsymbol{B}\coloneqq\bigg{(}\tilde{w}_{m}\left(\boldsymbol{x}_{j}-\boldsymbol{M}_{\boldsymbol{\sigma}}^{-1}\odot\boldsymbol{\ell}\right)\bigg{)}_{j=1,\,\l\in\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}}^{N}\
\in\mathbb{R}^{N\times|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|},$
(2.15)

where by definition (2.12) each row of $\boldsymbol{B}$ contains at most
$(2m+1)^{d}$ nonzeros. In doing so, the NFFT in Algorithm 2.1 can be
formulated in matrix-vector notation such that we receive the approximation
$\boldsymbol{A}\approx\boldsymbol{B}\boldsymbol{F}\boldsymbol{D}$ of (2.7),
cf. [53, p. 383]. In other words, the NFFT uses the approximation

$\displaystyle\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}}\approx\frac{1}{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|\cdot\hat{w}(\boldsymbol{k})}\sum_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}},m}(\boldsymbol{x}_{j})}\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}(\boldsymbol{M}_{\boldsymbol{\sigma}}^{-1}\odot\,\boldsymbol{\ell})}\,\tilde{w}_{m}\left(\boldsymbol{x}_{j}-\boldsymbol{M}_{\boldsymbol{\sigma}}^{-1}\odot\boldsymbol{\ell}\right).$

######  Remark 2.3.

It has to be pointed out that because of consistency the factor
$|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|^{-1}$ is here not
located in the matrix $\boldsymbol{F}$ as usual but in the matrix
$\boldsymbol{D}$.

###  2.2 The adjoint NFFT

Now we proceed with the adjoint problem (2.6). As already seen, this can be
written as $\boldsymbol{h}=\boldsymbol{A}^{*}\boldsymbol{f}$ with the adjoint
matrix $\boldsymbol{A}^{*}$ of (2.7). Thus, using the matrices (2.13), (2.14)
and (2.15) we receive the approximation
$\boldsymbol{A}^{*}\approx\boldsymbol{D}^{*}\boldsymbol{F}^{*}\boldsymbol{B}^{*}$,
such that a fast algorithm for the adjoint problem can be denoted as follows.

######  Algorithm 2.4 (adjoint NFFT).

For $d,N\in\mathbb{N}$ let $\boldsymbol{x}_{j}\in\mathbb{T}^{d}$,
$j=1,\dots,N$, be given points as well as $f_{j}\in\mathbb{C}$ given
coefficients. Further, we are given the oversampling factor $\sigma\geq 1$,
$2\mathbb{N}\ni M_{\sigma}\coloneqq 2\lceil\lceil\sigma M\rceil/2\rceil$,
$\boldsymbol{M}_{\boldsymbol{\sigma}}\coloneqq
M_{\sigma}\cdot\boldsymbol{1}_{d}$, as well as the window function $w$, the
truncated function $w_{m}$ with truncation parameter $m\ll M_{\sigma},$ and
their 1-periodic versions $\tilde{w}$ and $\tilde{w}_{m}$.

  1. 1.

Compute the sparse sums

$g_{\boldsymbol{\ell}}\coloneqq\sum_{j=1}^{N}f_{j}\,\tilde{w}_{m}\left(\boldsymbol{x}_{j}-\boldsymbol{M}_{\boldsymbol{\sigma}}^{-1}\odot\boldsymbol{\ell}\right),\quad\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}.$

$\mathcal{O}(N)$

  2. 2.

Compute

$\hat{g}_{\boldsymbol{k}}\coloneqq\frac{1}{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|}\sum_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}}g_{\boldsymbol{\ell}}\,\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{k}(\boldsymbol{M}_{\boldsymbol{\sigma}}^{-1}\odot\,\boldsymbol{\ell})},\quad\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}},$

by means of a $d$-variate FFT.
$\mathcal{O}(|\mathcal{I}_{\boldsymbol{M}}|\log(|\mathcal{I}_{\boldsymbol{M}}|))$

  3. 3.

Set

$\tilde{h}_{\boldsymbol{k}}\coloneqq\frac{\hat{g}_{\boldsymbol{k}}}{\hat{w}(\boldsymbol{k})},\quad\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}.$

$\mathcal{O}(|\mathcal{I}_{\boldsymbol{M}}|)$

Output: $\tilde{h}_{\boldsymbol{k}}\approx h_{\boldsymbol{k}}$,
$\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}$, cf. (2.6). Complexity:
$\mathcal{O}(|\mathcal{I}_{\boldsymbol{M}}|\log(|\mathcal{I}_{\boldsymbol{M}}|)+N)$

The algorithms presented in this section (Algorithms 2.1 and 2.4) are part of
the NFFT software [35]. For algorithmic details we refer to [36].

##  3 Direct inversion using density compensation

Having introduced the fast methods for nonequispaced data, we remark that
various applications such as MRI, solution of PDEs, etc. are interested in the
inverse problem, i. e., instead of the evaluation of (2.5) the aim is
computing the Fourier coefficients $\hat{f}_{\boldsymbol{k}}$,
$\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}$, from given nonequispaced data
$f(\boldsymbol{x}_{j})$, $j=1,\dots,N$. Therefore, this section shall be
dedicated to this task.

To clarify the major dissimilarity between equispaced and nonequispaced data,
we start considering the equispaced case. When evaluating at the points
$x_{\boldsymbol{j}}=\tfrac{1}{n}\boldsymbol{j}\in\mathbb{T}^{d}$,
$\boldsymbol{j}\in\mathcal{I}_{\boldsymbol{n}}$, with $\boldsymbol{n}\coloneqq
n\cdot\boldsymbol{1}_{d}$ and $|\mathcal{I}_{\boldsymbol{n}}|=N$, the
nonequispaced Fourier matrix
$\boldsymbol{A}\in\mathbb{C}^{N\times|\mathcal{I}_{\boldsymbol{M}}|}$ in (2.7)
turns into the equispaced Fourier matrix
$\boldsymbol{F}\in\mathbb{C}^{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|\times|\mathcal{I}_{\boldsymbol{M}}|}$
from (2.14) with $|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|=N$.
Thereby, it results from the geometric sum formula that

$\boldsymbol{F}^{*}\boldsymbol{F}=\Bigg{(}\sum_{\boldsymbol{j}\in\mathcal{I}_{\boldsymbol{n}}}\mathrm{e}^{2\pi\mathrm{i}(\boldsymbol{k}-\boldsymbol{\ell})\boldsymbol{j}/{n}}\Bigg{)}_{\boldsymbol{k},\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}}}=N\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|},\quad\text{
if }\,|\mathcal{I}_{\boldsymbol{M}}|\leq N,$ (3.1)

as well as

$\boldsymbol{F}\boldsymbol{F}^{*}=\Bigg{(}\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}(\boldsymbol{j}-\boldsymbol{h})/{n}}\Bigg{)}_{\boldsymbol{j},\boldsymbol{h}\in\mathcal{I}_{\boldsymbol{n}}}=|\mathcal{I}_{\boldsymbol{M}}|\cdot\boldsymbol{I}_{N},\quad\text{
if }\,|\mathcal{I}_{\boldsymbol{M}}|\geq N$ (3.2)

and $|\mathcal{I}_{\boldsymbol{M}}|$ is divisible by $N$. Thus, in the
equispaced setting a one-sided inverse is given by the (scaled) adjoint
matrix. However, when considering arbitrary points
$\boldsymbol{x}_{j}\in\mathbb{T}^{d}$, $j=1,\dots,N$, this property is lost,
i. e., for the nonequispaced Fourier matrix we have

$\displaystyle\boldsymbol{A}^{*}\boldsymbol{A}\neq
N\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}\quad\text{ and
}\quad\boldsymbol{A}\boldsymbol{A}^{*}\neq|\mathcal{I}_{\boldsymbol{M}}|\cdot\boldsymbol{I}_{N}.$
(3.3)

Because of this, more effort is needed in the nonequispaced setting.

In general, we face the following two problems.

  1. (1)

Solve the linear system

$\begin{split}\boldsymbol{A}\boldsymbol{\hat{f}}&=\boldsymbol{f},\end{split}$
(3.4)

i. e., reconstruct the Fourier coefficients
$\boldsymbol{\hat{f}}=(\hat{f}_{\boldsymbol{k}})_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}$
from given function values $\boldsymbol{f}=(f(\boldsymbol{x}_{j}))_{j=1}^{N}$.
This problem is referred to as inverse NDFT (iNDFT) and an efficient solver
shall be called inverse NFFT (iNFFT).

  2. (2)

Solve the linear system

$\begin{split}\boldsymbol{A}^{*}\boldsymbol{f}&=\boldsymbol{h},\end{split}$
(3.5)

i. e., reconstruct the coefficients $\boldsymbol{f}=(f_{j})_{j=1}^{N}$ from
given data
$\boldsymbol{h}=(h_{\boldsymbol{k}})_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}$.
This problem is referred to as inverse adjoint NDFT (iNDFT*) and an efficient
solver shall be called inverse adjoint NFFT (iNFFT*).

Note that in both problems the numbers $|\mathcal{I}_{\boldsymbol{M}}|$ and
$N$ are independent, such that the nonequispaced Fourier matrix
$\boldsymbol{A}\in\mathbb{C}^{N\times|\mathcal{I}_{\boldsymbol{M}}|}$ in (2.7)
is generally rectangular.

At first, we restrict our attention to problem (3.4). When considering
iterative inversion procedures as those mentioned in the introduction, these
methods require multiple iteration steps by definition. Therefore, multiple
matrix vector multiplications with the system matrix $\boldsymbol{A}$, or
rather multiple applications of the NFFT (see Algorithm 2.1), are needed to
compute a solution. To reduce the computational effort, we now proceed, in
contrast to this iterated procedures, with so-called direct methods. In the
setting of problem (3.4) we hereby mean methods, where for a fixed set of
points $\boldsymbol{x}_{j}$, $j=1,\dots,N$, the reconstruction of
$\boldsymbol{\hat{f}}$ from given $\boldsymbol{f}$ can be realized with the
same number of arithmetic operations as a single application of an adjoint
NFFT (see Algorithm 2.4). To achieve this, a certain precomputational step is
compulsory, since the adjoint NFFT does not yield an inversion of the NFFT per
se, see (3.3). Although this precomputations might be rather costly, they need
to be done only once for a given set of points $\boldsymbol{x}_{j}$,
$j=1,\dots,N$. In fact, the actual reconstruction step is very efficient.
Therefore, direct methods are especially beneficial in case we are given fixed
points for several measurement vectors $\boldsymbol{f}$.

In this section we focus on a direct inversion method for solving problem
(3.4) that utilizes so-called _sampling density compensation_. To this end, we
consider the integral (2.2) and introduce a corresponding quadrature formula.
In contrast to the already known equispaced approximation (2.3) we now assume
given arbitrary, nonequispaced points $\boldsymbol{x}_{j}\in\mathbb{T}^{d}$,
$j=1,\dots,N$. Thereby, the Fourier coefficients (2.2) are approximated by a
general quadrature rule using quadrature weights $w_{j}\in\mathbb{C}$,
$j=1,\dots,N$, which are needed for sampling density compensation due to the
nonequispaced sampling. Thus, for a trigonometric polynomial (2.8) we have

$\hat{f}_{\boldsymbol{k}}=c_{\boldsymbol{k}}(f)\approx
h_{\boldsymbol{k}}^{\mathrm{w}}\coloneqq\sum_{j=1}^{N}w_{j}\,f(\boldsymbol{x}_{j})\,\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}},\quad\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}.$
(3.6)

Using the nonequispaced Fourier matrix
$\boldsymbol{A}\in\mathbb{C}^{N\times|\mathcal{I}_{\boldsymbol{M}}|}$ in
(2.7), the diagonal matrix of weights
$\boldsymbol{W}\coloneqq\mathrm{diag}(w_{j})_{j=1}^{N}\in\mathbb{C}^{N\times
N}$ as well as the vector
$\boldsymbol{h}^{\mathrm{w}}\coloneqq(h_{\boldsymbol{k}}^{\mathrm{w}})_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}$,
the nonequispaced quadrature rule (3.6) can be written as
$\boldsymbol{\hat{f}}\approx\boldsymbol{h}^{\mathrm{w}}\coloneqq\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{f}$.
For achieving a fast computation method we make use of the approximation of
the adjoint NFFT, cf. Section 2.2, i. e., the final approximation is given by

$\boldsymbol{\hat{f}}\approx\boldsymbol{\tilde{h}}^{\mathrm{w}}\coloneqq\boldsymbol{D}^{*}\boldsymbol{F}^{*}\boldsymbol{B}^{*}\boldsymbol{W}\boldsymbol{f},$
(3.7)

with the matrices
$\boldsymbol{D}\in\mathbb{C}^{|\mathcal{I}_{\boldsymbol{M}}|\times|\mathcal{I}_{\boldsymbol{M}}|}$,
$\boldsymbol{F}\in\mathbb{C}^{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|\times|\mathcal{I}_{\boldsymbol{M}}|}$
and
$\boldsymbol{B}\in\mathbb{R}^{N\times|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|}$
defined in (2.13), (2.14) and (2.15). In other words, for density compensation
methods the already mentioned precomputations consist of computing the
quadrature weights $w_{j}\in\mathbb{C}$, $j=1,\dots,N$, while the actual
reconstruction step includes only one adjoint NFFT (see Algorithm 2.4) applied
to the scaled measurement vector $\boldsymbol{W}\boldsymbol{f}$.

The aim of all density compensation techniques is then to choose appropriate
weights $w_{j}\in\mathbb{C}$, $j=1,\dots,N$, such that the underlying
quadrature (3.6) is preferably exact. In the following we have a look at the
specific choice of the so-called _density compensation factors_ $w_{j}$.

An intuitive approach for density compensation is based on geometry, where
each sample is considered as representative of a certain surrounding area, as
in numerical integration. The weights for each sample can be obtained for
instance by constructing a Voronoi diagram and calculating the area of each
cell, see e. g. [58]. This approach of _Voronoi weights_ is well-known and
widely used in practice. However, it does not necessarily yield a good
approximation (3.7), which is why we examine some more sophisticated
approaches in the remainder of this section.

To this end, this section is organized as follows. Firstly, in Section 3.1 we
introduce density compensation factors $w_{j}$, $j=1,\dots,N,$ that lead to an
exact reconstruction formula (3.6) for all trigonometric polynomials (2.8) of
degree $\boldsymbol{M}$. In addition to the theoretical results, we also
discuss methods for the numerical computation. Secondly, in Section 3.2 we
show that it is reasonable to consider the inversion problem (3.4) and density
compensation via (3.7) for bandlimited functions $f\in
L_{1}(\mathbb{R}^{d})\cap C_{0}(\mathbb{R}^{d})$ as well. Subsequently, we
summarize our previous findings by presenting a general error bound on density
compensation factors in Section 3.3. Finally, in Section 3.4 we reconsider
certain approaches from literature and illustrate their connection among each
other as well as to the method introduced in Section 3.1.

######  Remark 3.1.

Recapitulating, we have a closer look at some possible interpretation
perspectives on the reconstruction (3.7).

  1. (i)

If we define $\boldsymbol{g}\coloneqq\boldsymbol{W}\boldsymbol{f}$, i. e.,
each entry of $\boldsymbol{f}$ is scaled with respect to the points
$\boldsymbol{x}_{j}$, $j=1,\dots,N$, the approximation (3.7) can be written as
$\boldsymbol{\hat{f}}\approx\boldsymbol{D}^{*}\boldsymbol{F}^{*}\boldsymbol{B}^{*}\boldsymbol{g}$.
As mentioned before, this coincides with an ordinary adjoint NFFT applied to a
modified coefficient vector $\boldsymbol{g}$.

  2. (ii)

By defining the matrix
$\boldsymbol{\tilde{B}}\coloneqq\boldsymbol{W}^{*}\boldsymbol{B}$, i. e.,
scaling the rows of $\boldsymbol{B}$ with respect to the points
$\boldsymbol{x}_{j}$, $j=1,\dots,N$, the approximation (3.7) can be written as
$\boldsymbol{\hat{f}}\approx\boldsymbol{D}^{*}\boldsymbol{F}^{*}\boldsymbol{\tilde{B}}^{*}\boldsymbol{f}$.
In this sense, density compensation can also be seen as a modification of the
adjoint NFFT and its application to the original coefficient vector.

Note that (i) is the common viewpoint. However, we keep (ii) in mind, since
this allows treating density compensation methods as an optimization of the
sparse matrix
$\boldsymbol{B}\in\mathbb{R}^{N\times|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|}$
in (2.15), as it shall be done in Section 4. We remark that density
compensation methods allow only $N$ degrees of freedom.

###  3.1 Exact quadrature weights for trigonometric polynomials

Similar to [28], we aim to introduce density compensation factors $w_{j}$,
$j=1,\dots,N,$ that lead to an exact reconstruction formula (3.6) for all
trigonometric polynomials (2.8) of degree $\boldsymbol{M}$. To this end, we
firstly examine certain properties that arise from (3.6) being exact.

######  Theorem 3.2.

Let a polynomial degree $\boldsymbol{M}\in(2\mathbb{N})^{d}$, nonequispaced
points $\boldsymbol{x}_{j}\in\mathbb{T}^{d}$, $j=1,\dots,N,$ and quadrature
weights $w_{j}\in\mathbb{C}$ be given. Then an exact reconstruction formula
(3.6) for trigonometric polynomials (2.8) with maximum degree $\boldsymbol{M}$
satisfying

$\displaystyle\hat{f}_{\boldsymbol{k}}=c_{\boldsymbol{k}}(f)=h_{\boldsymbol{k}}^{\mathrm{w}},\quad\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}},$
(3.8)

implies the following equivalent statements.

  1. (i)

The quadrature rule

$\int_{\mathbb{T}^{d}}f(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}=\sum_{j=1}^{N}w_{j}f(\boldsymbol{x}_{j})$
(3.9)

is exact for all trigonometric polynomials (2.8) with maximum degree
$\boldsymbol{M}$.

  2. (ii)

The linear system of equations

$\displaystyle\left[\boldsymbol{A}^{T}\boldsymbol{w}\right]_{\boldsymbol{k}}=\sum_{j=1}^{N}w_{j}\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}}=\delta_{\boldsymbol{0},\boldsymbol{k}}\,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}=\left\\{\begin{array}[]{ll}1&\colon\boldsymbol{k}=\boldsymbol{0}\\\
0&\colon\text{otherwise}\end{array}\right\\}},\quad\boldsymbol{k}\in{\mathcal{I}_{\boldsymbol{M}}},$
(3.12)

is fulfilled with the matrix
$\boldsymbol{A}\in\mathbb{C}^{N\times|\mathcal{I}_{\boldsymbol{M}}|}$ in (2.7)
and $\boldsymbol{w}\coloneqq\left(w_{j}\right)_{j=1}^{N}$.

###### Proof.

text

(3.8) $\Rightarrow$ (i): By inserting the definition (2.8) of a trigonometric
polynomial of degree $\boldsymbol{M}$ into the integral considered in (3.9) we
have

$\displaystyle\int_{\mathbb{T}^{d}}f(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}$
$\displaystyle=\sum_{\boldsymbol{k}\in{\mathcal{I}_{\boldsymbol{M}}}}\hat{f}_{\boldsymbol{k}}\cdot\int_{\mathbb{T}^{d}}\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}}\,\mathrm{d}\boldsymbol{x}=\sum_{\boldsymbol{k}\in{\mathcal{I}_{\boldsymbol{M}}}}\hat{f}_{\boldsymbol{k}}\cdot\delta_{\boldsymbol{0},\boldsymbol{k}}=\hat{f}_{\boldsymbol{0}},$
(3.13)

with the Kronecker delta $\delta_{\boldsymbol{0},\boldsymbol{k}}$. Now using
the property (3.8) as well as definition (3.6) of
$h_{\boldsymbol{k}}^{\mathrm{w}}$ we proceed with

$\displaystyle\hat{f}_{\boldsymbol{0}}=h_{\boldsymbol{0}}^{\mathrm{w}}$
$\displaystyle=\sum_{j=1}^{N}w_{j}f(\boldsymbol{x}_{j})\sum_{\boldsymbol{k}\in{\mathcal{I}_{\boldsymbol{M}}}}\mathrm{e}^{0}=\sum_{j=1}^{N}w_{j}f(\boldsymbol{x}_{j}),$
(3.14)

such that (3.13) combined with (3.14) yields the assertion (3.9).

(i) $\Rightarrow$ (ii): Inserting the definition (2.8) of a trigonometric
polynomial of degree $\boldsymbol{M}$ into the right-hand side of (3.9)
implies

$\displaystyle\sum_{j=1}^{N}w_{j}f(\boldsymbol{x}_{j})$
$\displaystyle=\sum_{j=1}^{N}w_{j}\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\hat{f}_{\boldsymbol{k}}\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}}=\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\hat{f}_{\boldsymbol{k}}\,\sum_{j=1}^{N}w_{j}\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}}.$
(3.15)

This together with the property (i) and (3.13) leads to

$\displaystyle\hat{f}_{\boldsymbol{0}}=\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\hat{f}_{\boldsymbol{k}}\,\sum_{j=1}^{N}w_{j}\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}}$

and thus to assertion (3.12).

(ii) $\Rightarrow$ (i): Combining (3.13), (3.12) and (3.15) yields the
assertion via

$\displaystyle\int_{\mathbb{T}^{d}}f(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}$
$\displaystyle=\sum_{\boldsymbol{k}\in{\mathcal{I}_{\boldsymbol{M}}}}\hat{f}_{\boldsymbol{k}}\cdot\delta_{\boldsymbol{0},\boldsymbol{k}}=\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\hat{f}_{\boldsymbol{k}}\,\sum_{j=1}^{N}w_{j}\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}}=\sum_{j=1}^{N}w_{j}f(\boldsymbol{x}_{j}).$

∎

######  Remark 3.3.

Comparable results can also be found in literature. A fundamental theorem in
numerical integration, see [70], states that for any integral
$\int_{\mathbb{T}^{d}}f(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}$ there
exists an exact quadrature rule (3.9), i. e., optimal points
$\boldsymbol{x}_{j}\in\mathbb{T}^{d}$ and weights $w_{j}\in\mathbb{C}$,
$j=1,\dots,N$, such that (3.9) is fulfilled. In [31, Lemma 2.6] it was shown
that for given points $\boldsymbol{x}_{j}\in\mathbb{T}^{d}$, $j=1,\dots,N,$
certain quadrature weights $w_{j}$ can be stated by means of frame theoretical
considerations which lead to an exact quadrature rule (3.9) by definition.
Moreover, it was shown (cf. [31, Lemma 3.6]) that these weights are the ones
with minimal (weighted) $\ell_{2}$-norm, which are already known under the
name “least squares quadrature”, see [34]. According to [34, Sec. 2.1] these
quadrature weights $w_{j}$, $j=1,\dots,N,$ can be found by solving a linear
system of equations $\boldsymbol{\Phi}\boldsymbol{w}=\boldsymbol{v}$, where
$\Phi_{\boldsymbol{k},j}=\phi_{\boldsymbol{k}}(\boldsymbol{x}_{j})$ and
$v_{\boldsymbol{k}}=\int_{\mathbb{T}^{d}}\phi_{\boldsymbol{k}}(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}$
for a given set of basis functions
$\\{\phi_{\boldsymbol{k}}\\}_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}$.
In our setting we have
$\phi_{\boldsymbol{k}}(\boldsymbol{x})=\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}}$
and therefore

$\displaystyle\boldsymbol{\Phi}=\left(\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}}\right)_{\boldsymbol{k},j}=\boldsymbol{A}^{T}\quad\text{
and }\quad
v_{\boldsymbol{k}}=\int_{\mathbb{T}^{d}}1\cdot\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}}\,\mathrm{d}\boldsymbol{x}=\delta_{\boldsymbol{0},\boldsymbol{k}},$

i. e., the same linear system of equations as in (3.12). We remark that both
[31] and [34] state the results in the case $d=1$, a generalization to $d>1$,
however, is straight-forward.

By means of Theorem 3.2 we can now give a condition that guaranties (3.6)
being exact for all trigonometric polynomials (2.8) with maximum degree
$\boldsymbol{M}$.

######  Corollary 3.4.

The two statements (i) and (ii) in Theorem 3.2 are not equivalent to property
(3.8), since (3.12) does not imply an exact reconstruction in (3.6).

However, an augmented variant of (3.12), namely

$\displaystyle\sum_{j=1}^{N}w_{j}\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}}=\delta_{\boldsymbol{0},\boldsymbol{k}},\quad\boldsymbol{k}\in{\mathcal{I}_{\boldsymbol{2M}}},$
(3.16)

yields an exact reconstruction
$\hat{f}_{\boldsymbol{k}}=h_{\boldsymbol{k}}^{\mathrm{w}}$ in (3.6) for
trigonometric polynomials (2.8) with maximum degree $\boldsymbol{M}$.
Additionally, (3.16) implies the matrix equation
$\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A}=\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}$
with $\boldsymbol{A}\in\mathbb{C}^{N\times|\mathcal{I}_{\boldsymbol{M}}|}$ in
(2.7) and the identity matrix
$\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}$ of size
$|\mathcal{I}_{\boldsymbol{M}}|$.

###### Proof.

Utilizing definitions (3.6) and (2.8) we have

$\displaystyle h_{\boldsymbol{k}}^{\mathrm{w}}$
$\displaystyle=\sum_{j=1}^{N}w_{j}\,\Bigg{(}\sum_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}}}\hat{f}_{\boldsymbol{\ell}}\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{\ell}\boldsymbol{x}_{j}}\Bigg{)}\,\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}}=\sum_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}}}\hat{f}_{\boldsymbol{\ell}}\,\sum_{j=1}^{N}w_{j}\,\mathrm{e}^{2\pi\mathrm{i}(\boldsymbol{\ell}-\boldsymbol{k})\boldsymbol{x}_{j}}$
$\displaystyle=\sum_{\begin{subarray}{c}\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}}\\\
(\boldsymbol{\ell}-\boldsymbol{k})\in\mathcal{I}_{\boldsymbol{M}}\end{subarray}}\hat{f}_{\boldsymbol{\ell}}\,\sum_{j=1}^{N}w_{j}\,\mathrm{e}^{2\pi\mathrm{i}(\boldsymbol{\ell}-\boldsymbol{k})\boldsymbol{x}_{j}}+\sum_{\begin{subarray}{c}\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}}\\\
(\boldsymbol{\ell}-\boldsymbol{k})\notin\mathcal{I}_{\boldsymbol{M}}\end{subarray}}\hat{f}_{\boldsymbol{\ell}}\,\sum_{j=1}^{N}w_{j}\,\mathrm{e}^{2\pi\mathrm{i}(\boldsymbol{\ell}-\boldsymbol{k})\boldsymbol{x}_{j}},\quad\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}.$

Since (3.12) only holds for
$\boldsymbol{k},\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}}$ with
$(\boldsymbol{\ell}-\boldsymbol{k})\in\mathcal{I}_{\boldsymbol{M}}$, this
implies

$\displaystyle h_{\boldsymbol{k}}^{\mathrm{w}}$
$\displaystyle=\hat{f}_{\boldsymbol{k}}+\sum_{\begin{subarray}{c}\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}}\\\
(\boldsymbol{\ell}-\boldsymbol{k})\notin\mathcal{I}_{\boldsymbol{M}}\end{subarray}}\hat{f}_{\boldsymbol{\ell}}\,\sum_{j=1}^{N}w_{j}\,\mathrm{e}^{2\pi\mathrm{i}(\boldsymbol{\ell}-\boldsymbol{k})\boldsymbol{x}_{j}},\quad\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}},$

where for all
$\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}\setminus\\{\boldsymbol{0}\\}$
there exists an $\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}}$ with
$(\boldsymbol{\ell}-\boldsymbol{k})\in\mathcal{I}_{\boldsymbol{2M}}\setminus\mathcal{I}_{\boldsymbol{M}}$.

As $(\boldsymbol{\ell}-\boldsymbol{k})\in\mathcal{I}_{\boldsymbol{2M}}$ for
$\boldsymbol{k},\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}}$, the
augmented variant (3.16) yields

$\displaystyle h_{\boldsymbol{k}}^{\mathrm{w}}$
$\displaystyle=\sum_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}}}\hat{f}_{\boldsymbol{\ell}}\,\sum_{j=1}^{N}w_{j}\,\mathrm{e}^{2\pi\mathrm{i}(\boldsymbol{\ell}-\boldsymbol{k})\boldsymbol{x}_{j}}=\sum_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}}}\hat{f}_{\boldsymbol{k}}\cdot\delta_{\boldsymbol{0},\boldsymbol{k}}=\hat{f}_{\boldsymbol{k}},\quad\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}.$

Moreover, since
$\delta_{(\boldsymbol{\ell}-\boldsymbol{k}),\boldsymbol{0}}=\delta_{\boldsymbol{k},\boldsymbol{\ell}}$,
the condition (3.16) implies

$\displaystyle\delta_{\boldsymbol{k},\boldsymbol{\ell}}=\sum_{j=1}^{N}w_{j}\,\mathrm{e}^{2\pi\mathrm{i}(\boldsymbol{\ell}-\boldsymbol{k})\boldsymbol{x}_{j}}=\sum_{j=1}^{N}\,\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}}\left(w_{j}\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{\ell}\boldsymbol{x}_{j}}\right),\quad\boldsymbol{k},\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}}.$

In matrix-vector notation this can be written as
$\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A}=\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}$
with $\boldsymbol{A}\in\mathbb{C}^{N\times|\mathcal{I}_{\boldsymbol{M}}|}$ in
(2.7) and the identity matrix
$\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}$ of size
$|\mathcal{I}_{\boldsymbol{M}}|$. We remark that this matrix equation
immediately shows that we have an exact reconstruction of the form (3.8),
since if
$\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A}=\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}$
is fulfilled, (3.4) implies that
$\boldsymbol{\hat{f}}=\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A}\boldsymbol{\hat{f}}=\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{f}$.
∎

######  Remark 3.5.

Let $f\in L_{2}(\mathbb{T}^{d})$ be an arbitrary 1-periodic function (2.1).
Then (3.16) yields

$\displaystyle h_{\boldsymbol{k}}^{\mathrm{w}}$
$\displaystyle=\sum_{\boldsymbol{\ell}\in\mathbb{Z}^{d}}c_{\boldsymbol{\ell}}(f)\,\sum_{j=1}^{N}w_{j}\,\mathrm{e}^{2\pi\mathrm{i}(\boldsymbol{\ell}-\boldsymbol{k})\boldsymbol{x}_{j}}$
$\displaystyle=\sum_{\begin{subarray}{c}\boldsymbol{\ell}\in\mathbb{Z}^{d}\\\
(\boldsymbol{\ell}-\boldsymbol{k})\in\mathcal{I}_{\boldsymbol{2M}}\end{subarray}}c_{\boldsymbol{\ell}}(f)\,\sum_{j=1}^{N}w_{j}\,\mathrm{e}^{2\pi\mathrm{i}(\boldsymbol{\ell}-\boldsymbol{k})\boldsymbol{x}_{j}}+\sum_{\begin{subarray}{c}\boldsymbol{\ell}\in\mathbb{Z}^{d}\\\
(\boldsymbol{\ell}-\boldsymbol{k})\notin\mathcal{I}_{\boldsymbol{2M}}\end{subarray}}c_{\boldsymbol{\ell}}(f)\,\sum_{j=1}^{N}w_{j}\,\mathrm{e}^{2\pi\mathrm{i}(\boldsymbol{\ell}-\boldsymbol{k})\boldsymbol{x}_{j}}$
$\displaystyle=c_{\boldsymbol{k}}(f)+\sum_{\begin{subarray}{c}\boldsymbol{\ell}\in\mathbb{Z}^{d}\\\
(\boldsymbol{\ell}-\boldsymbol{k})\notin\mathcal{I}_{\boldsymbol{2M}}\end{subarray}}c_{\boldsymbol{\ell}}(f)\,\sum_{j=1}^{N}w_{j}\,\mathrm{e}^{2\pi\mathrm{i}(\boldsymbol{\ell}-\boldsymbol{k})\boldsymbol{x}_{j}},\quad\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}},$

i. e., for a function $f\in L_{2}(\mathbb{T}^{d})$ we only have a good
approximation in case the coefficients $c_{\boldsymbol{\ell}}(f)$ are small
for $\boldsymbol{\ell}\notin\mathcal{I}_{\boldsymbol{M}}$, whereas this
reconstruction can only be exact for $f$ being a trigonometric polynomial
(2.8).

####  3.1.1 Practical computation in the underdetermined setting
$|\mathcal{I}_{\boldsymbol{2M}}|\leq N$

So far, we have seen in Corollary 3.4 that an exact solution
$\boldsymbol{w}=(w_{j})_{j=1}^{N}$ to the linear system (3.16) leads to an
exact reconstruction formula (3.6) for all trigonometric polynomials (2.8)
with maximum degree $\boldsymbol{M}$. Therefore, we aim to use this condition
(3.16) to numerically find optimal density compensation factors
$w_{j}\in\mathbb{C}$, $j=1,\dots,N$.

Having a closer look at the condition (3.16) we recognize that it can be
written as the linear system of equations
$\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}^{T}\,\boldsymbol{w}=\boldsymbol{e}_{\boldsymbol{0}}$
with the matrix
$\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}\in\mathbb{C}^{N\times|\mathcal{I}_{\boldsymbol{2M}}|}$,
cf. (2.7), and right side
$\boldsymbol{e}_{\boldsymbol{0}}\coloneqq\left(\delta_{\boldsymbol{0},\boldsymbol{k}}\right)_{\boldsymbol{k}\in{\mathcal{I}_{\boldsymbol{2M}}}}$.
We remark that in contrast to
$\boldsymbol{A}\in\mathbb{C}^{N\times|\mathcal{I}_{\boldsymbol{M}}|}$ we now
deal with the enlarged matrix
$\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}\in\mathbb{C}^{N\times|\mathcal{I}_{\boldsymbol{2M}}|}$,
such that single matrix operations are more costly. Nevertheless, Corollary
3.4 yields a direct inversion method for (3.4), where the system
$\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}^{T}\,\boldsymbol{w}=\boldsymbol{e}_{\boldsymbol{0}}$
needs to be solved only once for fixed points
$\boldsymbol{x}_{j}\in\mathbb{T}^{d}$, $j=1,\dots,N$. Its solution
$\boldsymbol{w}$ can then be used to efficiently approximate
$\boldsymbol{\hat{f}}$ for multiple measurement vectors $\boldsymbol{f}$,
whereas iterative methods for (3.4) need to solve
$\boldsymbol{A}\boldsymbol{\hat{f}}=\boldsymbol{f}$ each time.

As already mentioned in [34, Sec. 3.1] an exact solution to (3.16) can only be
found if $|\mathcal{I}_{\boldsymbol{2M}}|\leq N$, i. e., in case
$\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}^{T}\,\boldsymbol{w}=\boldsymbol{e}_{\boldsymbol{0}}$
is an underdetermined system of equations. By [34, Lem. 3.1] this system has
at least one solution, which is why we may choose the one with minimal
$\ell_{2}$-norm. If
$\mathrm{rank}(\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|})={|\mathcal{I}_{\boldsymbol{2M}}|}$,
then the system
$\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}^{T}\,\boldsymbol{w}=\boldsymbol{e}_{\boldsymbol{0}}$
is consistent and the unique solution is given by the normal equations of
second kind

$\displaystyle\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}^{T}\overline{\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}}\,\boldsymbol{v}=\boldsymbol{e}_{\boldsymbol{0}},\quad\overline{\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}}\,\boldsymbol{v}=\boldsymbol{w}.$
(3.17)

More precisely, we may compute the vector $\boldsymbol{v}$ using an iterative
procedure such as the CG algorithm, such that only matrix multiplications with
$\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}^{T}$ and
$\overline{\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}}$ are needed.
Since fast multiplication with
$\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}^{T}$ and
$\overline{\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}}$ can easily be
realized by means of an adjoint NFFT (see Algorithm 2.4) and an NFFT (see
Algorithm 2.1), respectively, computing the solution $\boldsymbol{w}$ to
(3.17) is of complexity
$\mathcal{O}(|\mathcal{I}_{\boldsymbol{2M}}|\log(|\mathcal{I}_{\boldsymbol{2M}}|)+N)$,
where

$|\mathcal{I}_{\boldsymbol{2M}}|=(2M)^{d}=2^{d}M^{d}=2^{d}\,|\mathcal{I}_{\boldsymbol{M}}|.$

Thus, in order to receive exact quadrature weights $w_{j}\in\mathbb{C}$,
$j=1,\dots,N,$ via (3.17) we need to satisfy the full rank condition
$\mathrm{rank}(\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|})={|\mathcal{I}_{\boldsymbol{2M}}|}$.
In case of a low rank matrix
$\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}\in\mathbb{C}^{N\times|\mathcal{I}_{\boldsymbol{2M}}|}$
for $|\mathcal{I}_{\boldsymbol{2M}}|\leq N$, we may still use (3.17) to obtain
a least squares approximation to (3.16).

####  3.1.2 Practical computation in the overdetermined setting
$|\mathcal{I}_{\boldsymbol{2M}}|>N$

In the setting $|\mathcal{I}_{\boldsymbol{2M}}|>N$, we cannot expect finding
an exact solution $\boldsymbol{w}$ to (3.16), since we have to deal with an
overdetermined system possessing more conditions than variables. However, we
still aim to numerically find optimal density compensation factors
$w_{j}\in\mathbb{C}$, $j=1,\dots,N$, by considering a least squares
approximation to (3.16) that minimizes
$\big{\|}\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}^{T}\,\boldsymbol{w}-\boldsymbol{e}_{\boldsymbol{0}}\big{\|}_{2}$.
In [8, Thm. 1.1.2] it was shown that every least squares solution satisfies
the normal equations of first kind

$\displaystyle\overline{\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}}\,\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}^{T}\,\boldsymbol{w}=\overline{\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}}\,\boldsymbol{e}_{\boldsymbol{0}}.$
(3.18)

By means of definitions of
$\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}\in\mathbb{C}^{N\times|\mathcal{I}_{\boldsymbol{2M}}|}$,
cf. (2.7), and
$\boldsymbol{e}_{\boldsymbol{0}}=\left(\delta_{\boldsymbol{0},\boldsymbol{k}}\right)_{\boldsymbol{k}\in{\mathcal{I}_{\boldsymbol{2M}}}}$
we simplify the right hand side via

$\displaystyle\overline{\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}}\,\boldsymbol{e}_{\boldsymbol{0}}=\Bigg{(}\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{2M}}}\delta_{\boldsymbol{0},\boldsymbol{k}}\,\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}}\Bigg{)}_{j=1}^{N}=\boldsymbol{1}_{N}.$

Since fast multiplication with
$\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}^{T}$ and
$\overline{\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}}$ can easily be
realized by means of an adjoint NFFT (see Algorithm 2.4) and an NFFT (see
Algorithm 2.1), respectively, the solution $\boldsymbol{w}$ to (3.18) can be
computed iteratively by means of the CG algorithm in
$\mathcal{O}(|\mathcal{I}_{\boldsymbol{2M}}|\log(|\mathcal{I}_{\boldsymbol{2M}}|)+N)$
arithmetic operations. Note that the solution to (3.18) is only unique if the
full rank condition
$\mathrm{rank}(\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|})=N$ is
satisfied, cf. [8, p. 7]. We remark that the computed weight matrix
$\boldsymbol{W}=\operatorname*{diag}(\boldsymbol{w})$ can further be used in
an iterative procedure as in [53, Alg. 7.27] to improve the approximation of
$\boldsymbol{\hat{f}}$.

The previous considerations lead to the following algorithms.

######  Algorithm 3.6 (Computation of the optimal density compensation
factors).

For $d,N\in\mathbb{N}$ let $\boldsymbol{x}_{j}\in\mathbb{T}^{d}$,
$j=1,\dots,N$, as well as $M\in 2\mathbb{N}$ and $\boldsymbol{M}\coloneqq
M\cdot\boldsymbol{1}_{d}$ be given.

  1. 1.

Compute $|\mathcal{I}_{\boldsymbol{2M}}|=(2M)^{d}$. $\mathcal{O}(1)$

  2. 2.

If $|\mathcal{I}_{\boldsymbol{2M}}|\leq N$

     * Compute the solution $\boldsymbol{v}$ to (3.17) iteratively using the NFFT.
.
$\mathcal{O}(|\mathcal{I}_{\boldsymbol{2M}}|\log(|\mathcal{I}_{\boldsymbol{2M}}|)+N)$

     * Compute the solution $\boldsymbol{w}=\overline{\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}}\,\boldsymbol{v}$, see (3.17), using an NFFT.
.
$\mathcal{O}(|\mathcal{I}_{\boldsymbol{2M}}|\log(|\mathcal{I}_{\boldsymbol{2M}}|)+N)$

elseif $|\mathcal{I}_{\boldsymbol{2M}}|>N$

     * Compute the solution $\boldsymbol{w}$ to (3.18) iteratively using the NFFT.
.
$\mathcal{O}(|\mathcal{I}_{\boldsymbol{2M}}|\log(|\mathcal{I}_{\boldsymbol{2M}}|)+N)$

  3. 3.

Compose
$\boldsymbol{W}=\operatorname*{diag}(\boldsymbol{w})\in\mathbb{C}^{N\times
N}$. $\mathcal{O}(N)$

Output: weights matrix $\boldsymbol{W}$ Complexity:
$\mathcal{O}(|\mathcal{I}_{\boldsymbol{2M}}|\log(|\mathcal{I}_{\boldsymbol{2M}}|)+N)$

######  Algorithm 3.7 (iNFFT – density compensation approach).

For $d,N\in\mathbb{N}$ let $\boldsymbol{x}_{j}\in\mathbb{T}^{d}$,
$j=1,\dots,N$, as well as $\boldsymbol{f}\in\mathbb{C}^{N}$, $M\in
2\mathbb{N}$ and $\boldsymbol{M}\coloneqq M\cdot\boldsymbol{1}_{d}$ be given.

  1. 0.

Precompute the weights matrix $\boldsymbol{W}$ using Algorithm 3.6.

  2. 1.

Compute
$\boldsymbol{\tilde{h}}^{\mathrm{w}}\coloneqq\boldsymbol{D}^{*}\boldsymbol{F}^{*}\boldsymbol{B}^{*}\boldsymbol{W}\boldsymbol{f}$,
cf. (3.7), by means of an adjoint NFFT.

Output:
$\boldsymbol{\tilde{h}}^{\mathrm{w}}\approx\boldsymbol{\hat{f}}\in\mathbb{C}^{|\mathcal{I}_{\boldsymbol{M}}|}$,
cf. (3.4). Complexity:
$\mathcal{O}(|\mathcal{I}_{\boldsymbol{M}}|\log(|\mathcal{I}_{\boldsymbol{M}}|)+N)$

###  3.2 Bandlimited functions

In some numerical examples, such as in MRI, we are concerned with bandlimited
functions $f\in L_{1}(\mathbb{R}^{d})\cap C_{0}(\mathbb{R}^{d})$ instead of
trigonometric polynomials $f\in L_{2}(\mathbb{T}^{d})$ in (2.8), cf. [15]. In
the following we show that it is reasonable to consider the inversion problem
(3.4) as well as density compensation via (3.7) for bandlimited functions as
well.

To this end, let $f\in L_{1}(\mathbb{R}^{d})\cap C_{0}(\mathbb{R}^{d})$ be a
bandlimited function with bandwidth $\boldsymbol{M}$, i. e., its (continuous)
Fourier transform

$\hat{f}(\boldsymbol{v})\coloneqq\int\limits_{\mathbb{R}^{d}}f(\boldsymbol{x})\,\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{v}\boldsymbol{x}}\,\mathrm{d}\boldsymbol{x},\quad\boldsymbol{v}\in\mathbb{R}^{d},$
(3.19)

is supported on $\left[-\tfrac{M}{2},\tfrac{M}{2}\right)^{d}$. Utilizing this
fact, we have $\hat{f}\in L_{1}(\mathbb{R}^{d})$ and thus by the Fourier
inversion theorem [53, Thm. 2.10] the inverse Fourier transform of $f$ can be
written as

$f(\boldsymbol{x})=\int\limits_{\mathbb{R}^{d}}\hat{f}(\boldsymbol{v})\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{v}\boldsymbol{x}}\,\mathrm{d}\boldsymbol{v}=\int\limits_{\left[-\frac{M}{2},\frac{M}{2}\right)^{d}}\hat{f}(\boldsymbol{v})\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{v}\boldsymbol{x}}\,\mathrm{d}\boldsymbol{v},\quad\boldsymbol{x}\in\mathbb{R}^{d}.$
(3.20)

Analogous to (2.3), the approximation using equispaced quadrature points
$\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}$ yields

$f(\boldsymbol{x})\approx\frac{\left(\frac{M}{2}-(-\frac{M}{2})\right)^{d}}{|\mathcal{I}_{\boldsymbol{M}}|}\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\hat{f}(\boldsymbol{k})\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}}=\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\hat{f}(\boldsymbol{k})\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}},\quad\boldsymbol{x}\in\mathbb{R}^{d},$
(3.21)

such that evaluation at the given nonequispaced points
$\boldsymbol{x}_{j}\in\left[-\tfrac{1}{2},\tfrac{1}{2}\right)^{d}$,
$j=1,\dots,N$, leads to

$f(\boldsymbol{x}_{j})\approx\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\hat{f}(\boldsymbol{k})\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}}.$

By means of the definition (2.7) of the matrix
$\boldsymbol{A}\in\mathbb{C}^{N\times|\mathcal{I}_{\boldsymbol{M}}|}$ this can
be written as $\boldsymbol{f}\approx\boldsymbol{A}\boldsymbol{\hat{f}}$, where
we used the notation
$\boldsymbol{\hat{f}}\coloneqq(\hat{f}(\boldsymbol{k}))_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}$
in this setting. Thus, also for bandlimited functions $f$ its evaluations at
points $\boldsymbol{x}_{j}$ can be approximated in the form (2.5), such that
it is reasonable to consider the inversion problem (3.4) for bandlimited
functions as well.

Considering (3.19) we are given an exact formula for the evaluation of the
Fourier transform $\hat{f}$. However, in practical applications, such as MRI,
this is only a hypothetical case, since $f$ cannot be sampled on whole
$\mathbb{R}^{d}$, cf. [15]. Due to a limited coverage of space by the
acquisition, the function $f$ is typically only known on a bounded domain,
w.l.o.g. for $\boldsymbol{x}\in\left[-\tfrac{1}{2},\tfrac{1}{2}\right)^{d}$.
Thus, we have to deal with the approximation

$\hat{f}(\boldsymbol{v})\approx\int\limits_{\left[-\tfrac{1}{2},\tfrac{1}{2}\right)^{d}}f(\boldsymbol{x})\,\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{v}\boldsymbol{x}}\,\mathrm{d}\boldsymbol{x},\quad\boldsymbol{v}\in\left[-\tfrac{M}{2},\tfrac{M}{2}\right)^{d}.$
(3.22)

Using the nonequispaced quadrature rule in (3.6), we find that evaluation at
uniform grid points $\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}$ can be
approximated via

$\hat{f}(\boldsymbol{k})\approx\tilde{h}(\boldsymbol{k})\coloneqq\sum_{j=1}^{N}w_{j}\,f(\boldsymbol{x}_{j})\,\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}},\quad\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}.$

This is to say, equispaced samples of the Fourier transform of a bandlimited
function may be approximated in the same form
$(\hat{f}(\boldsymbol{k}))_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}=\boldsymbol{\hat{f}}\approx\boldsymbol{h}^{\mathrm{w}}\coloneqq\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{f}$
as in (3.6), where we used the notation
$\boldsymbol{h}^{\mathrm{w}}\coloneqq(\tilde{h}(\boldsymbol{k}))_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}$
in this setting. Moreover, we extend this approximation onto the whole
interval $\left[-\tfrac{M}{2},\tfrac{M}{2}\right)^{d}$, i. e., we consider

$\hat{f}(\boldsymbol{v})\approx\tilde{h}(\boldsymbol{v})\coloneqq\sum_{j=1}^{N}w_{j}\,f(\boldsymbol{x}_{j})\,\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{v}\boldsymbol{x}_{j}},\quad\boldsymbol{v}\in\left[-\tfrac{M}{2},\tfrac{M}{2}\right)^{d}.$
(3.23)

So all in all, we have seen that it is reasonable to study the inversion
problem (3.4) and the associated density compensation via (3.7) for
bandlimited functions as well. Analogous to Section 3.1 we now aim to find a
numerical method for computing suitable weights $w_{j}\in\mathbb{C}$,
$j=1,\dots,N$, such that the reconstruction formula (3.23) is preferably
exact. To this end, we have a closer look at (3.23) being exact and start with
analogous considerations as in Theorem 3.2.

######  Theorem 3.8.

Let a bandwidth $\boldsymbol{M}\in\mathbb{N}^{d}$, nonequispaced points
$\boldsymbol{x}_{j}\in\left[-\tfrac{1}{2},\tfrac{1}{2}\right)^{d}$ as well as
quadrature weights $w_{j}\in\mathbb{C}$, $j=1,\dots,N,$ be given. Then an
exact reconstruction formula (3.23) for bandlimited functions $f\in
L_{1}(\mathbb{R}^{d})\cap C_{0}(\mathbb{R}^{d})$ with bandwidth
$\boldsymbol{M}$, i. e.,

$\hat{f}(\boldsymbol{v})=\tilde{h}(\boldsymbol{v})=\sum_{j=1}^{N}w_{j}f(\boldsymbol{x}_{j})\,\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{v}\boldsymbol{x}_{j}},\quad\boldsymbol{v}\in\left[-\tfrac{M}{2},\tfrac{M}{2}\right)^{d},$
(3.24)

implies that the quadrature rule

$\displaystyle\int_{\mathbb{R}^{d}}f(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}=\sum_{j=1}^{N}w_{j}f(\boldsymbol{x}_{j})$

is exact for all bandlimited functions $f\in L_{1}(\mathbb{R}^{d})\cap
C_{0}(\mathbb{R}^{d})$ with bandwidth $\boldsymbol{M}$.

###### Proof.

By (3.19) the assumption (3.24) can be written as

$\int\limits_{\mathbb{R}^{d}}f(\boldsymbol{x})\,\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{v}\boldsymbol{x}}\,\mathrm{d}\boldsymbol{x}=\hat{f}(\boldsymbol{v})=\sum_{j=1}^{N}w_{j}f(\boldsymbol{x}_{j})\,\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{v}\boldsymbol{x}_{j}},\quad\boldsymbol{v}\in\left[-\tfrac{M}{2},\tfrac{M}{2}\right)^{d}.$
(3.25)

Especially, for $\boldsymbol{v}=\boldsymbol{0}$ evaluation of (3.25) yields
the assertion

$\int\limits_{\mathbb{R}^{d}}f(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}=\hat{f}(\boldsymbol{0})=\sum_{j=1}^{N}w_{j}f(\boldsymbol{x}_{j}).$

∎

However, in contrast to Theorem 3.2, this Theorem 3.8 does not yield an
explicit condition for computing suitable weights $w_{j}\in\mathbb{C}$,
$j=1,\dots,N$. To derive a numerical procedure anyway, we generalize the
notion of an exact reconstruction $\tilde{h}$ of $f$ and have a look at the
theory of tempered distributions. To this end, let
$\mathscr{S}(\mathbb{R}^{d})$ be the Schwartz space of rapidly decaying
functions, cf. [53, Sec. 4.2.1]. The tempered Dirac distribution $\delta$
shall be defined by
$\langle\delta,\varphi\rangle\coloneqq\int_{\mathbb{R}^{d}}\varphi(\boldsymbol{v})\,\delta(\boldsymbol{v})\,\mathrm{d}\boldsymbol{v}=\varphi(\boldsymbol{0})$
for all $\varphi\in\mathscr{S}(\mathbb{R}^{d})$, cf. [53, Ex. 4.36]. For a
slowly increasing function $f\colon\mathbb{R}^{d}\to\mathbb{C}$ satisfying
$|f(\boldsymbol{x})|\leq c(1+\|\boldsymbol{x}\|_{2})^{n}$ almost everywhere
with $c>0$ and $n\in\mathbb{N}_{0}$, the induced distribution $T_{f}$ shall be
defined by $\langle
T_{f},\varphi\rangle\coloneqq\int_{\mathbb{R}^{d}}\varphi(\boldsymbol{x})\,f(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}$
for all $\varphi\in\mathscr{S}(\mathbb{R}^{d})$. For a detailed introduction
to the topic we refer to [53, Sec. 4.2.1 and Sec. 4.3].

Then the following property can be shown.

######  Theorem 3.9.

Let nonequispaced points
$\boldsymbol{x}_{j}\in\left[-\frac{1}{2},\frac{1}{2}\right)^{d}$,
$j=1,\dots,N,$ and quadrature weights $w_{j}\in\mathbb{C}$ be given. Further
let $T_{f}$ be the distribution induced by some bandlimited function $f\in
L_{1}(\mathbb{R}^{d})\cap C_{0}(\mathbb{R}^{d})$ with bandwidth
$\boldsymbol{M}$. Then

$\displaystyle\langle\delta,\varphi\rangle=\langle
T_{\xi},\varphi\rangle,\quad\varphi\in\mathscr{S}(\mathbb{R}^{d}),$ (3.26)

with

$\displaystyle\xi(\boldsymbol{v})\coloneqq\sum_{j=1}^{N}w_{j}\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{v}\boldsymbol{x}_{j}},\quad\boldsymbol{v}\in\mathbb{R}^{d},$
(3.27)

implies

$\displaystyle\langle\hat{T}_{f},\varphi\rangle=\langle
T_{\tilde{h}},\varphi\rangle,\quad\varphi\in\mathscr{S}(\mathbb{R}^{d}),$
(3.28)

with the function $\tilde{h}$ defined in (3.23).

###### Proof.

Using the definition of the function $\tilde{h}$ in (3.23) as well as the fact
that the inversion formula (3.20) holds for all
$\boldsymbol{x}\in\mathbb{R}^{d}$, we have

$\displaystyle\langle T_{\tilde{h}},\varphi\rangle$
$\displaystyle=\int_{\mathbb{R}^{d}}\varphi(\boldsymbol{v})\,\sum_{j=1}^{N}w_{j}f(\boldsymbol{x}_{j})\,\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{v}\boldsymbol{x}_{j}}\,\mathrm{d}\boldsymbol{v}$
$\displaystyle=\int_{\mathbb{R}^{d}}\varphi(\boldsymbol{v})\,\sum_{j=1}^{N}w_{j}\left(\int_{\mathbb{R}^{d}}\hat{f}(\boldsymbol{u})\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{u}\boldsymbol{x}_{j}}\,\mathrm{d}\boldsymbol{u}\right)\,\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{v}\boldsymbol{x}_{j}}\,\mathrm{d}\boldsymbol{v}$
$\displaystyle=-\int_{\mathbb{R}^{d}}\hat{f}(\boldsymbol{u})\int_{\mathbb{R}^{d}}\varphi(\boldsymbol{u}-\boldsymbol{v})\,\sum_{j=1}^{N}w_{j}\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{v}\boldsymbol{x}_{j}}\,\mathrm{d}\boldsymbol{v}\,\mathrm{d}\boldsymbol{u}.$

Hence, by (3.26) this implies

$\displaystyle\langle T_{\tilde{h}},\varphi\rangle$
$\displaystyle=-\int_{\mathbb{R}^{d}}\hat{f}(\boldsymbol{u})\int_{\mathbb{R}^{d}}\varphi(\boldsymbol{u}-\boldsymbol{v})\,\delta(\boldsymbol{v})\,\mathrm{d}\boldsymbol{v}\,\mathrm{d}\boldsymbol{u}$
$\displaystyle=\int_{\mathbb{R}^{d}}\hat{f}(\boldsymbol{u})\int_{\mathbb{R}^{d}}\varphi(\boldsymbol{v})\,\delta(\boldsymbol{u}-\boldsymbol{v})\,\mathrm{d}\boldsymbol{v}\,\mathrm{d}\boldsymbol{u}=\int_{\mathbb{R}^{d}}\hat{f}(\boldsymbol{u})\,\varphi(\boldsymbol{u})\,\mathrm{d}\boldsymbol{u}.$

∎

Considering the property (3.28), we remark that this indeed states an exact
reconstruction $\hat{f}=\tilde{h}$ in the sense of tempered distributions, as
distinct from (3.24). Since it is known by Corollary 3.4 that the condition
(3.16) yields an exact reconstruction for trigonometric polynomials, we aim to
use this result to compute suitable weights $w_{j}\in\mathbb{C}$,
$j=1,\dots,N$, for bandlimited functions as well. To this end, suppose we have
(3.16), i. e.,

$\displaystyle\sum_{j=1}^{N}w_{j}\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}}=\delta_{\boldsymbol{0},\boldsymbol{k}},\quad\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{2M}}.$

Then this yields

$\displaystyle\varphi(\boldsymbol{0})=\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{2M}}}\varphi(\boldsymbol{k})\,\sum_{j=1}^{N}w_{j}\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}},\quad\varphi\in\mathscr{S}(\mathbb{R}^{d}).$
(3.29)

Having a look at Theorem 3.9, an exact reconstruction (3.28) is implied by
(3.26), i. e.,

$\displaystyle\varphi(\boldsymbol{0})=\langle\delta,\varphi\rangle=\langle
T_{\xi},\varphi\rangle=\int_{\mathbb{R}^{d}}\varphi(\boldsymbol{v})\,\sum_{j=1}^{N}w_{j}\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{v}\boldsymbol{x}_{j}}\,\mathrm{d}\boldsymbol{v},\quad\varphi\in\mathscr{S}(\mathbb{R}^{d}).$

Thus, the property (3.29) that is fulfilled by (3.16) could be interpreted as
an equispaced quadrature of (3.26) at integer frequencies
$\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{2M}}$.

######  Remark 3.10.

We remark that for deriving the quadrature rule (3.29) from (3.26), we
implicitly truncate the integral bounds in (3.26) as

$\displaystyle\langle
T_{\xi},\varphi\rangle=\int_{\mathbb{R}^{d}}\varphi(\boldsymbol{v})\,\sum_{j=1}^{N}w_{j}\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{v}\boldsymbol{x}_{j}}\,\mathrm{d}\boldsymbol{v}$
$\displaystyle\approx\int\limits_{[-M,M)^{d}}\varphi(\boldsymbol{v})\,\sum_{j=1}^{N}w_{j}\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{v}\boldsymbol{x}_{j}}\,\mathrm{d}\boldsymbol{v}$
$\displaystyle=\int_{\mathbb{R}^{d}}\varphi(\boldsymbol{v})\,\sum_{j=1}^{N}w_{j}\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{v}\boldsymbol{x}_{j}}\,\chi_{\left[-M,M\right)^{d}}(\boldsymbol{v})\,\mathrm{d}\boldsymbol{v},$

i. e., instead of (3.27) we rather deal with a distribution induced by

$\displaystyle\tilde{\xi}(\boldsymbol{v})\coloneqq\xi(\boldsymbol{v})\,\chi_{\left[-M,M\right)^{d}}(\boldsymbol{v}),\quad\boldsymbol{v}\in\mathbb{R}^{d}.$
(3.30)

However, an analogous implication as in Theorem 3.9 cannot be shown when using
the function (3.30) instead of (3.27).

As seen in Remark 3.10, our numerical method for computing suitable weights
$w_{j}\in\mathbb{C}$, $j=1,\dots,N$, can also be derived by means of a
quadrature formula applied to the property
$\langle\delta,\hat{\varphi}\rangle=\langle
T_{\tilde{\xi}},\hat{\varphi}\rangle$ with $\tilde{\xi}$ defined in (3.30).
Having a closer look at this property, the following equivalent
characterization can be shown.

######  Theorem 3.11.

Let a bandwidth $\boldsymbol{M}\in\mathbb{N}^{d}$, nonequispaced points
$\boldsymbol{x}_{j}\in\left[-\frac{1}{2},\frac{1}{2}\right)^{d}$ and
quadrature weights $w_{j}\in\mathbb{C}$, $j=1,\dots,N,$ be given. Then the
following two statements are equivalent.

  1. (ii)

For all $\varphi\in\mathscr{S}(\mathbb{R}^{d})$ we have
$\langle\delta,\hat{\varphi}\rangle=\langle
T_{\tilde{\xi}},\hat{\varphi}\rangle$ with $\tilde{\xi}$ defined in (3.30).

  2. (iii)

We have $\langle 1,\varphi\rangle=\langle T_{\psi},\varphi\rangle$ for all
$\varphi\in\mathscr{S}(\mathbb{R}^{d})$, where

$\displaystyle\psi(\boldsymbol{x})\coloneqq\sum_{j=1}^{N}w_{j}\cdot|\mathcal{I}_{\boldsymbol{2M}}|\,\mathrm{sinc}\left(2M\pi\left(\boldsymbol{x}_{j}-\boldsymbol{x}\right)\right),\quad\boldsymbol{x}\in\mathbb{R}^{d},$

with the $d$-variate $\mathrm{sinc}$ function
$\mathrm{sinc}(\boldsymbol{x})\coloneqq\prod_{t=1}^{d}\mathrm{sinc}(x_{t})$
and

$\mathrm{sinc}(x)\coloneqq\left\\{\begin{array}[]{ll}\frac{\sin x}{x}&\quad
x\in\mathbb{R}\setminus\\{0\\}\,,\\\\[4.30554pt] 1&\quad
x=0\,.\end{array}\right.$

###### Proof.

By the definition $\langle\hat{T},\varphi\rangle=\langle
T,\hat{\varphi}\rangle$ of the Fourier transform of a tempered distribution
$T\in\mathscr{S}^{\prime}(\mathbb{R}^{d})$ we have $\langle
1,\varphi\rangle=\langle\delta,\hat{\varphi}\rangle$, cf. [53, Ex. 4.46].
Moreover, the distribution induced by (3.30) can be rewritten using the
Fourier transform (3.19) as

$\displaystyle\langle T_{\tilde{\xi}},\hat{\varphi}\rangle$
$\displaystyle=\int_{\mathbb{R}^{d}}\hat{\varphi}(\boldsymbol{v})\,\tilde{\xi}(\boldsymbol{v})\,\mathrm{d}\boldsymbol{v}=\sum_{j=1}^{N}w_{j}\,\int\limits_{\left[-M,M\right)^{d}}\hat{\varphi}(\boldsymbol{v})\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{v}\boldsymbol{x}_{j}}\,\mathrm{d}\boldsymbol{v}$
$\displaystyle=\sum_{j=1}^{N}w_{j}\int\limits_{\left[-M,M\right)^{d}}\left(\int_{\mathbb{R}^{d}}\varphi(\boldsymbol{x})\,\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{v}\boldsymbol{x}}\,\mathrm{d}\boldsymbol{x}\right)\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{v}\boldsymbol{x}_{j}}\,\mathrm{d}\boldsymbol{v}$
$\displaystyle=\sum_{j=1}^{N}w_{j}\int_{\mathbb{R}^{d}}\varphi(\boldsymbol{x})\int\limits_{\left[-M,M\right)^{d}}\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{v}(\boldsymbol{x}_{j}-\boldsymbol{x})}\,\mathrm{d}\boldsymbol{v}\,\mathrm{d}\boldsymbol{x}.$
(3.31)

The inner integral can be determined by

$\displaystyle\int_{\left[-M,M\right)^{d}}\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{v}(\boldsymbol{x}_{j}-\boldsymbol{x})}\,\mathrm{d}\boldsymbol{v}=|\mathcal{I}_{\boldsymbol{2M}}|\,\mathrm{sinc}\left(2M\pi\left(\boldsymbol{x}_{j}-\boldsymbol{x}\right)\right),$

such that (3.2) shows the equality $\langle
T_{\tilde{\xi}},\hat{\varphi}\rangle=\langle T_{\psi},\varphi\rangle$. Hence,
the assertions (ii) and (iii) are equivalent. ∎

Thus, since the statements (ii) and (iii) of Theorem 3.11 are equivalent, one
could also consider the property (iii) for deriving a numerical method to
compute suitable weights $w_{j}\in\mathbb{C}$, $j=1,\dots,N$. To this end, we
have a closer look at

$\displaystyle\int_{\mathbb{R}^{d}}\varphi(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}=\langle
1,\varphi\rangle=\langle
T_{\psi},\varphi\rangle=\int_{\mathbb{R}^{d}}\varphi(\boldsymbol{x})\,\sum_{j=1}^{N}w_{j}\,\cdot|\mathcal{I}_{\boldsymbol{2M}}|\,\mathrm{sinc}\left(2M\pi\left(\boldsymbol{x}_{j}-\boldsymbol{x}\right)\right)\,\mathrm{d}\boldsymbol{x}$
(3.32)

for all $\varphi\in\mathscr{S}(\mathbb{R}^{d})$. Due to the integrals on both
sides of (3.32), we need to discretize twice and therefore use the same
quadrature rule on both sides of (3.32). For better comparability to (3.29) we
utilize the same number $|\mathcal{I}_{\boldsymbol{2M}}|$ of equispaced
quadrature points
$\boldsymbol{y}_{\boldsymbol{\ell}}\coloneqq(\boldsymbol{2M})^{-1}\odot\boldsymbol{\ell}$,
$\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{2M}}$, as in (3.29), i. e., we
consider

$\displaystyle\sum_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{2M}}}\varphi(\boldsymbol{y}_{\boldsymbol{\ell}})=\sum_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{2M}}}\varphi(\boldsymbol{y}_{\boldsymbol{\ell}})\,\sum_{j=1}^{N}w_{j}\cdot|\mathcal{I}_{\boldsymbol{2M}}|\,\mathrm{sinc}\left(2M\pi\left(\boldsymbol{x}_{j}-\boldsymbol{y}_{\boldsymbol{\ell}}\right)\right).$

In order that this applies for all $\varphi\in\mathscr{S}(\mathbb{R}^{d})$, we
need to satisfy

$\displaystyle
1=\sum_{j=1}^{N}w_{j}\cdot|\mathcal{I}_{\boldsymbol{2M}}|\,\mathrm{sinc}\left(2M\pi\left(\boldsymbol{x}_{j}-\boldsymbol{y}_{\boldsymbol{\ell}}\right)\right),\quad\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{2M}},$
(3.33)

i. e., one could also compute weights $w_{j}\in\mathbb{C}$, $j=1,\dots,N$, as
a least squares solution to the linear system of equations (3.33). Hence, it
merely remains the comparison of the two computation schemes.

######  Remark 3.12.

Since we derived discretizations out of both statements of Theorem 3.11, we
examine if also the two linear systems (3.16) and (3.33) are related.
Considering the statements in Theorem 3.11 we notice that in some sense they
are the Fourier transformed versions of each other. To this end, we need to
Fourier transform one of the linear systems for better comparability. More
precisely, we apply an iFFT of length $|\mathcal{I}_{\boldsymbol{2M}}|$, cf.
(2.4), to both sides of equation (3.33). Since the left side transforms to

$\displaystyle\sum_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{2M}}}1\cdot\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{y}_{\boldsymbol{\ell}}}=|\mathcal{I}_{\boldsymbol{2M}}|\cdot\delta_{\boldsymbol{0},\boldsymbol{k}},\quad\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{2M}},$

we obtain the transformed system

$\displaystyle\delta_{\boldsymbol{0},\boldsymbol{k}}=\sum_{j=1}^{N}w_{j}\sum_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{2M}}}\,\mathrm{sinc}\left(2M\pi\left(\boldsymbol{x}_{j}-\boldsymbol{y}_{\boldsymbol{\ell}}\right)\right)\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{y}_{\boldsymbol{\ell}}},\quad\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{2M}}.$
(3.34)

Comparing this linear system of equations to (3.16), we recognize an identical
structure. Hence, we have a closer look at the connection between the
expressions $\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}}$ and
$\sum_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{2M}}}\mathrm{sinc}\left(2M\pi\left(\boldsymbol{x}_{j}-\boldsymbol{y}_{\boldsymbol{\ell}}\right)\right)\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{y}_{\boldsymbol{\ell}}}$.

For this purpose, we consider the function
$f(\boldsymbol{t})=\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{t}\boldsymbol{x}}$,
$\boldsymbol{t}\in\left[-M,M\right)^{d}$, for fixed
$\boldsymbol{x}\in\mathbb{C}^{d}$. By means of
$\tilde{f}(\boldsymbol{t})\coloneqq\sum_{\boldsymbol{k}\in\mathbb{Z}^{d}}f(\boldsymbol{t}+{2M}\boldsymbol{k})$
we extend it into a $({2M})$-periodic function. This periodized version then
possesses the Fourier coefficients

$\displaystyle c_{\boldsymbol{\ell}}(\tilde{f})$
$\displaystyle=\frac{1}{|\mathcal{I}_{\boldsymbol{2M}}|}\int\limits_{\left[-M,M\right)^{d}}f(\boldsymbol{t})\,\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{t}\boldsymbol{y}_{\boldsymbol{\ell}}}\,\mathrm{d}\boldsymbol{t}$
$\displaystyle=\frac{1}{|\mathcal{I}_{\boldsymbol{2M}}|}\int\limits_{\left[-M,M\right)^{d}}\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{t}(\boldsymbol{x}-\boldsymbol{y}_{\boldsymbol{\ell}})}\,\mathrm{d}\boldsymbol{t}=\mathrm{sinc}\left(2M\pi\left(\boldsymbol{x}-\boldsymbol{y}_{\boldsymbol{\ell}}\right)\right),\quad\boldsymbol{\ell}\in\mathbb{Z}^{d},$

cf. (2.2), i. e., the Fourier expansion of  $\tilde{f}(\boldsymbol{t})$,
$\boldsymbol{t}\in\left[-M,M\right)^{d}$, for fixed $\boldsymbol{x}$ is given
by

$\displaystyle\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{t}\boldsymbol{x}}=\sum_{\boldsymbol{\ell}\in\mathbb{Z}^{d}}\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{t}\boldsymbol{y}_{\boldsymbol{\ell}}}\,\mathrm{sinc}\left(2M\pi\left(\boldsymbol{x}-\boldsymbol{y}_{\boldsymbol{\ell}}\right)\right),\quad\boldsymbol{x}\in\mathbb{C}^{d},$

cf. (2.1). Since  $\tilde{f}(\boldsymbol{t})$ is continuous and piecewise
differentiable, this Fourier series converges absolutely and uniformly, cf.
[49,  Ex. 1.22]. Thereby, we may consider the point evaluations at
$\boldsymbol{x}=\boldsymbol{x}_{j}$, $j=1,\dots,N$, and
$\boldsymbol{t}=\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{2M}}$, such that we
obtain the representation

$\displaystyle\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}}=\sum_{\boldsymbol{\ell}\in\mathbb{Z}^{d}}\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{y}_{\boldsymbol{\ell}}}\,\mathrm{sinc}\left(2M\pi\left(\boldsymbol{x}_{j}-\boldsymbol{y}_{\boldsymbol{\ell}}\right)\right).$
(3.35)

Thus, we recognize that (3.34) is a truncated version of (3.35). In other
words, this implies that the linear system (3.16) is equivalent to a
discretization of (3.32) incorporating infinitely many points
$\boldsymbol{y}_{\boldsymbol{\ell}}\in\mathbb{R}^{d}$ in (3.33).

######  Remark 3.13.

For bandlimited functions several fast evaluation methods including the
$\mathrm{sinc}$ function are known. The classical sampling theorem of Shannon-
Whittaker-Kotelnikov, see [72, 64, 42], states that any bandlimited function
$f\in L_{1}(\mathbb{R}^{d})\cap C_{0}(\mathbb{R}^{d})$ with maximum bandwidth
$\boldsymbol{M}$ can be recovered from its uniform samples
$f(\boldsymbol{L}^{-1}\odot{\boldsymbol{\ell}})$,
$\boldsymbol{\ell}\in\mathbb{Z}^{d}$, with $L\geq M$, $\boldsymbol{L}\coloneqq
L\cdot\boldsymbol{1}_{d}$, and we have

$\displaystyle
f(\boldsymbol{x})=\sum_{\boldsymbol{\ell}\in\mathbb{Z}^{d}}f(\boldsymbol{L}^{-1}\odot{\boldsymbol{\ell}})\,\mathrm{sinc}\left(L\pi\left(\boldsymbol{x}-\boldsymbol{L}^{-1}\odot{\boldsymbol{\ell}}\right)\right),\quad\boldsymbol{x}\in\mathbb{R}^{d}.$
(3.36)

Since the practical use of this sampling theorem is limited due to the
infinite number of samples, which is impossible in practice, and the very slow
decay of the sinc function, various authors such as [57, 68, 50, 47, 40]
considered the regularized Shannon sampling formula with localized sampling

$\displaystyle
f(\boldsymbol{x})\approx\sum_{\boldsymbol{\ell}\in\mathbb{Z}^{d}}f(\boldsymbol{L}^{-1}\odot{\boldsymbol{\ell}})\,\mathrm{sinc}\left(L\pi\left(\boldsymbol{x}-\boldsymbol{L}^{-1}\odot{\boldsymbol{\ell}}\right)\right)\,\varphi_{m}\left(\boldsymbol{x}-\boldsymbol{L}^{-1}\odot{\boldsymbol{\ell}}\right),\quad\boldsymbol{x}\in\mathbb{R}^{d},$
(3.37)

instead. Here $\varphi_{m}\colon\mathbb{R}^{d}\to[0,1]$ is a compactly
supported window function with truncation parameter
$m\in\mathbb{N}\setminus\\{1\\}$, such that for $\varphi_{m}$ with small
support the direct evaluation of (3.37) is efficient, see [40] for the
relation to the NFFT window functions.

On the other hand, in the one-dimensional setting a fast $\mathrm{sinc}$
transform was introduced in [39], which is based on the Clenshaw-Curtis
quadrature

$\displaystyle\mathrm{sinc}\left(L\pi\left(x-\tfrac{\ell}{L}\right)\right)\approx\sum_{k=0}^{n}w_{k}\,\mathrm{e}^{-\pi\mathrm{i}L(x-\frac{\ell}{L})z_{k}}$

using Chebyshev points $z_{k}=\cos(\frac{k\pi}{n})\in[-1,1]$, $k=0,\dots,n$,
and corresponding Clenshaw-Curtis weights $w_{k}>0$. Thereby, sums of the form

$\displaystyle
h(x)=\sum_{\ell\in\mathcal{I}_{T}}f\big{(}\tfrac{\ell}{L}\big{)}\,\mathrm{sinc}\left(L\pi\left(x-\tfrac{\ell}{L}\right)\right)\approx\sum_{k=0}^{n}w_{k}\,\bigg{(}\sum_{\ell\in\mathcal{I}_{T}}f\big{(}\tfrac{\ell}{L}\big{)}\,\mathrm{e}^{\pi\mathrm{i}\ell
z_{k}}\bigg{)}\,\mathrm{e}^{-\pi\mathrm{i}Lxz_{k}}$

with uniform truncation parameter $T\in 2\mathbb{N}$, can efficiently be
approximated by means of fast Fourier transforms. More precisely, for the term
in brackets one may utilize an NFFT, cf. (2.5). Then the resulting outer sum
can be computed using an NNFFT, also referred to as NFFT of type III, see [17,
20, 46] or [53, pp. 394–397].

###  3.3 General error bound

In this section we summarize our previous findings by presenting a general
error bound on density compensation factors computed by means of (3.16), that
applies to trigonometric polynomials, 1-periodic functions $f\in
L_{2}(\mathbb{T}^{d})\cap C(\mathbb{T}^{d})$ and bandlimited functions $f\in
L_{1}(\mathbb{R}^{d})\cap C_{0}(\mathbb{R}^{d})$ as well.

######  Theorem 3.14.

Let $p,q\in\\{1,2,\infty\\}$ with $\frac{1}{p}+\frac{1}{q}=1$. For given
$d,N\in\mathbb{N}$, $\boldsymbol{M}\in(2\mathbb{N})^{d}$ and nonequispaced
points $\boldsymbol{x}_{j}\in\left[-\tfrac{1}{2},\tfrac{1}{2}\right)^{d}$,
$j=1,\dots,N$, let
$\boldsymbol{A}\in\mathbb{C}^{N\times|\mathcal{I}_{\boldsymbol{M}}|}$ be the
nonequispaced Fourier matrix in (2.7). Further assume we can compute density
compensation factors
$\boldsymbol{W}=\operatorname*{diag}\left(w_{j}\right)_{j=1}^{N}\in\mathbb{C}^{N\times
N}$ by means of Algorithm 3.6, such that

$\displaystyle\sum_{j=1}^{N}w_{j}\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}}=\delta_{\boldsymbol{0},\boldsymbol{k}}+\varepsilon_{\boldsymbol{k}},\quad\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{2M}},$
(3.38)

with small $\varepsilon_{\boldsymbol{k}}\in\mathbb{R}$ for all
$\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{2M}}$.

Then there exists an $\varepsilon\geq 0$ such that the corresponding density
compensation procedure with $\boldsymbol{W}=\mathrm{diag}(w_{j})_{j=1}^{N}$
satisfies the following error bounds.

  1. (i)

For any trigonometric polynomial $f\in L_{2}(\mathbb{T}^{d})$ of degree
$\boldsymbol{M}$ given in (2.8) we have

$\displaystyle\big{\|}\boldsymbol{\hat{f}}-\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{f}\big{\|}_{p}$
$\displaystyle\leq|\mathcal{I}_{\boldsymbol{M}}|\,\varepsilon\cdot\big{\|}\boldsymbol{\hat{f}}\big{\|}_{p},$
(3.39)

where
$\boldsymbol{\hat{f}}\coloneqq(\hat{f}_{\boldsymbol{k}})_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}$
are the coefficients given in (2.8).

  2. (ii)

For any 1-periodic function $f\in L_{2}(\mathbb{T}^{d})\cap C(\mathbb{T}^{d})$
we have

$\displaystyle\big{\|}\boldsymbol{\hat{f}}-\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{f}\big{\|}_{p}$
$\displaystyle\leq|\mathcal{I}_{\boldsymbol{M}}|\,\varepsilon\cdot\big{\|}\boldsymbol{\hat{f}}\big{\|}_{p}+(N\,|\mathcal{I}_{\boldsymbol{M}}|)^{1/p}\,\|\boldsymbol{w}\|_{q}\cdot\|f-p_{\boldsymbol{M}}\|_{C(\mathbb{T}^{d})},$
(3.40)

where
$\boldsymbol{\hat{f}}\coloneqq(c_{\boldsymbol{k}}(f))_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}$
are the first $|\mathcal{I}_{\boldsymbol{M}}|$ coefficients given in (2.1) and
$p_{\boldsymbol{M}}$ is the best approximating trigonometric polynomial of
degree $\boldsymbol{M}$ of $f$.

  3. (iii)

For any bandlimited function $f\in L_{1}(\mathbb{R}^{d})\cap
C_{0}(\mathbb{R}^{d})$ with bandwidth $\boldsymbol{M}$ we have

$\displaystyle\big{\|}\boldsymbol{\hat{f}}-\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{f}\big{\|}_{p}$
$\displaystyle\leq|\mathcal{I}_{\boldsymbol{M}}|\,\varepsilon\cdot\big{\|}\boldsymbol{\hat{f}}\big{\|}_{p}+(N\,|\mathcal{I}_{\boldsymbol{M}}|)^{1/p}\,\|\boldsymbol{w}\|_{q}\cdot\|Q\|_{C(\mathbb{T}^{d})},$
(3.41)

where
$\boldsymbol{\hat{f}}\coloneqq(\hat{f}(\boldsymbol{k}))_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}$
are the integer evaluations of (3.19) and $Q$ in (3.47) is the pointwise
quadrature error of the equispaced quadrature rule (3.21).

###### Proof.

We start with some general considerations that are independent of the function
$f$. By (3.38) we can find
$\varepsilon\coloneqq\max_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}|\varepsilon_{\boldsymbol{k}}|\geq
0$, such that $|\varepsilon_{\boldsymbol{k}}|\leq\varepsilon$,
$\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{2M}}$, and thereby

$\displaystyle\left|\,\sum_{j=1}^{N}w_{j}\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}}-\delta_{\boldsymbol{0},\boldsymbol{k}}\,\right|\leq\varepsilon,\quad\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{2M}}.$

Then for all $\boldsymbol{k},\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}}$
with $(\boldsymbol{\ell}-\boldsymbol{k})\in\mathcal{I}_{\boldsymbol{2M}}$ this
yields

$\displaystyle\left|\left[\boldsymbol{E}_{\mathrm{r}}\right]_{\boldsymbol{k},\boldsymbol{\ell}}\,\right|=\left|\,\sum_{j=1}^{N}w_{j}\,\mathrm{e}^{2\pi\mathrm{i}(\boldsymbol{\ell}-\boldsymbol{k})\boldsymbol{x}_{j}}-\delta_{\boldsymbol{k},\boldsymbol{\ell}}\,\right|\leq\varepsilon,$

where
$\boldsymbol{E}_{\mathrm{r}}\coloneqq\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A}-\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}$.
Hence, we have

$\displaystyle\big{\|}\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A}-\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}\big{\|}_{1}=\max_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}}}\,\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\left|\left[\boldsymbol{E}_{\mathrm{r}}\right]_{\boldsymbol{k},\boldsymbol{\ell}}\,\right|\leq\max_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}}}\,\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\varepsilon=|\mathcal{I}_{\boldsymbol{M}}|\,\varepsilon,$
(3.42)
$\displaystyle\big{\|}\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A}-\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}\big{\|}_{\infty}=\max_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\,\sum_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}}}\left|\left[\boldsymbol{E}_{\mathrm{r}}\right]_{\boldsymbol{k},\boldsymbol{\ell}}\,\right|\leq\max_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\,\sum_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}}}\varepsilon=|\mathcal{I}_{\boldsymbol{M}}|\,\varepsilon,$
(3.43)

and

$\displaystyle\big{\|}\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A}-\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}\big{\|}_{\mathrm{F}}=\sqrt{\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\,\sum_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}}}\left|\left[\boldsymbol{E}_{\mathrm{r}}\right]_{\boldsymbol{k},\boldsymbol{\ell}}\,\right|^{2}}\leq\sqrt{\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\,\sum_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}}}\varepsilon^{2}}=|\mathcal{I}_{\boldsymbol{M}}|\,\varepsilon.$
(3.44)

Considering the approximation error of (3.6), it can be estimated by

$\displaystyle\big{\|}\boldsymbol{\hat{f}}-\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{f}\big{\|}_{p}$
$\displaystyle\leq\big{\|}\boldsymbol{\hat{f}}-\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A}\boldsymbol{\hat{f}}\big{\|}_{p}+\big{\|}\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A}\boldsymbol{\hat{f}}-\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{f}\big{\|}_{p}$
$\displaystyle=\big{\|}\big{(}\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A}-\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}\big{)}\boldsymbol{\hat{f}}\big{\|}_{p}+\big{\|}\boldsymbol{A}^{*}\boldsymbol{W}\big{(}\boldsymbol{A}\boldsymbol{\hat{f}}-\boldsymbol{f}\big{)}\big{\|}_{p}$
$\displaystyle\leq\big{\|}\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A}-\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}\big{\|}_{p}\cdot\big{\|}\boldsymbol{\hat{f}}\big{\|}_{p}+\big{\|}\boldsymbol{A}^{*}\boldsymbol{W}\big{\|}_{p}\cdot\big{\|}\boldsymbol{A}\boldsymbol{\hat{f}}-\boldsymbol{f}\big{\|}_{p}.$
(3.45)

Using $\boldsymbol{A}\in\mathbb{C}^{N\times|\mathcal{I}_{\boldsymbol{M}}|}$
from (2.7) as well as
$\boldsymbol{W}=\mathrm{diag}(w_{j})_{j=1}^{N}=\operatorname*{diag}(\boldsymbol{w})$
we have

$\displaystyle\big{\|}\boldsymbol{A}^{*}\boldsymbol{W}\big{\|}_{1}=\max_{j=1,\dots,N}\,\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}|w_{j}|\cdot\big{|}\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}}\big{|}\leq\max_{j=1,\dots,N}|w_{j}|\cdot\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}1=|\mathcal{I}_{\boldsymbol{M}}|\cdot\|\boldsymbol{w}\|_{\infty},$
$\displaystyle\big{\|}\boldsymbol{A}^{*}\boldsymbol{W}\big{\|}_{\infty}=\max_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\,\sum_{j=1}^{N}|w_{j}|\cdot\big{|}\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}}\big{|}\leq\sum_{j=1}^{N}|w_{j}|\cdot\max_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}1=\|\boldsymbol{w}\|_{1},$

and

$\displaystyle\big{\|}\boldsymbol{A}^{*}\boldsymbol{W}\big{\|}_{\mathrm{F}}=\sqrt{\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\,\sum_{j=1}^{N}|w_{j}|^{2}\cdot\big{|}\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}}\big{|}^{2}}\leq\sqrt{\sum_{j=1}^{N}|w_{j}|^{2}\cdot|\mathcal{I}_{\boldsymbol{M}}|}=\sqrt{|\mathcal{I}_{\boldsymbol{M}}|}\cdot\|\boldsymbol{w}\|_{2}.$

Hence, from (3.42) – (3.45) and $\|\cdot\|_{2}\leq\|\cdot\|_{\mathrm{F}}$ it
follows that

$\displaystyle\big{\|}\boldsymbol{\hat{f}}-\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{f}\big{\|}_{p}$
$\displaystyle\leq|\mathcal{I}_{\boldsymbol{M}}|\,\varepsilon\cdot\big{\|}\boldsymbol{\hat{f}}\big{\|}_{p}+\big{\|}\boldsymbol{A}\boldsymbol{\hat{f}}-\boldsymbol{f}\big{\|}_{p}\cdot|\mathcal{I}_{\boldsymbol{M}}|^{1/p}\,\|\boldsymbol{w}\|_{q}$
(3.46)

for $p\in\\{1,2,\infty\\}$ with $\frac{1}{p}+\frac{1}{q}=1$. Now it merely
remains to estimate
$\big{\|}\boldsymbol{A}\boldsymbol{\hat{f}}-\boldsymbol{f}\big{\|}_{p}$ for
the specific choice of $f$.

(i): Since a trigonometric polynomial (2.8) of degree $\boldsymbol{M}$
satisfies $\boldsymbol{A}\boldsymbol{\hat{f}}=\boldsymbol{f}$, the second
error term in (3.46) vanishes and we obtain the assertion (3.39).

(ii): When considering a general 1-periodic function $f\in
L_{2}(\mathbb{T}^{d})\cap C(\mathbb{T}^{d})$ in (2.1) we have

$\displaystyle\Big{|}\big{[}\boldsymbol{A}\boldsymbol{\hat{f}}-\boldsymbol{f}\big{]}_{j}\Big{|}$
$\displaystyle=\Bigg{|}f(\boldsymbol{x}_{j})-\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}c_{\boldsymbol{k}}(f)\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}}\Bigg{|}$
$\displaystyle\leq\max_{\boldsymbol{x}\in\mathbb{T}^{d}}\Bigg{|}\sum_{\boldsymbol{k}\in\mathbb{Z}^{d}\setminus\mathcal{I}_{\boldsymbol{M}}}c_{\boldsymbol{k}}(f)\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}}\Bigg{|}=\|f-p_{\boldsymbol{M}}\|_{C(\mathbb{T}^{d})}\quad
j=1,\dots,N,$

with the best approximating trigonometric polynomial $p_{\boldsymbol{M}}$ of
degree $\boldsymbol{M}$ of $f$. Thus, this yields
$\big{\|}\boldsymbol{A}\boldsymbol{\hat{f}}-\boldsymbol{f}\big{\|}_{p}\leq
N^{1/p}\,\|f-p_{\boldsymbol{M}}\|_{C(\mathbb{T}^{d})}$ and by (3.46) the
assertion (3.40).

(iii): For a bandlimited function $f\in L_{1}(\mathbb{R}^{d})\cap
C_{0}(\mathbb{R}^{d})$ with bandwidth $\boldsymbol{M}$ we may use the notation
$\boldsymbol{\hat{f}}\coloneqq(\hat{f}(\boldsymbol{k}))_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}$
as well as the inverse Fourier transform (3.20) to estimate

$\displaystyle\Big{|}\big{[}\boldsymbol{A}\boldsymbol{\hat{f}}-\boldsymbol{f}\big{]}_{j}\Big{|}$
$\displaystyle=\Bigg{|}f(\boldsymbol{x}_{j})-\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\hat{f}(\boldsymbol{k})\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}}\Bigg{|}\leq\max_{\boldsymbol{x}\in\mathbb{T}^{d}}|Q(\boldsymbol{x})|=\|Q\|_{C(\mathbb{T}^{d})},\quad
j=1,\dots,N,$

with the pointwise quadrature error

$\displaystyle
Q(\boldsymbol{x})\coloneqq\int_{\left[-\frac{M}{2},\frac{M}{2}\right)^{d}}\hat{f}(\boldsymbol{v})\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{v}\boldsymbol{x}}\,\mathrm{d}\boldsymbol{v}-\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\hat{f}(\boldsymbol{k})\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}}$
(3.47)

of the uniform quadrature rule (3.21). For detailed investigations of
quadrature errors for bandlimited functions we refer to [39, 27]. Hence, we
obtain
$\big{\|}\boldsymbol{A}\boldsymbol{\hat{f}}-\boldsymbol{f}\big{\|}_{p}\leq
N^{1/p}\,\|Q\|_{C(\mathbb{T}^{d})}$ and by (3.46) the assertion (3.41).

∎

By Corollary 3.4 it is known that in the setting of trigonometric polynomials
there is a linkage between an exact reconstruction (3.8) and the matrix
product $\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A}$ being equal to
identity $\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}$. The following
theorem shows that the error of the reconstruction (3.6) also affects the
condition of the matrix $\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A}$.

######  Theorem 3.15.

Let $\boldsymbol{A}\in\mathbb{C}^{N\times|\mathcal{I}_{\boldsymbol{M}}|}$ from
(2.7), $\boldsymbol{W}=\mathrm{diag}(w_{j})_{j=1}^{N}$ and $\varepsilon\geq 0$
be given as in Theorem 3.14. If additionally
$\varepsilon\,|\mathcal{I}_{\boldsymbol{M}}|<1$ is fulfilled, then we have

$\displaystyle
1\leq\kappa_{2}(\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A})\leq\frac{1+\varepsilon\,|\mathcal{I}_{\boldsymbol{M}}|}{1-\varepsilon\,|\mathcal{I}_{\boldsymbol{M}}|}$
(3.48)

for the condition number
$\kappa_{2}(\boldsymbol{X})\coloneqq\|\boldsymbol{X}\|_{2}\|\boldsymbol{X}^{-1}\|_{2}$.

###### Proof.

To estimate the condition number
$\kappa_{2}(\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A})$ we need to
determine the norms
$\big{\|}\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A}\big{\|}_{2}$ and
$\big{\|}(\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A})^{-1}\big{\|}_{2}$.
By (3.38) it is known that
$\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A}=\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}+\boldsymbol{\mathcal{E}}$,
where
$\boldsymbol{\mathcal{E}}\coloneqq\left(\varepsilon_{\boldsymbol{\ell}-\boldsymbol{k}}\right)_{\boldsymbol{\ell},\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}$,
and therefore we have

$\displaystyle\big{\|}\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A}\big{\|}_{2}=\big{\|}\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}+\boldsymbol{\mathcal{E}}\big{\|}_{2}\leq\big{\|}\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}\big{\|}_{2}+\|\boldsymbol{\mathcal{E}}\|_{2}.$
(3.49)

Moreover, it is known by the theory of Neumann series, cf. [67, Thm. 4.20],
that if
$\big{\|}\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}-\boldsymbol{T}\big{\|}_{2}<1$
holds for a matrix
$\boldsymbol{T}\in\mathbb{C}^{|\mathcal{I}_{\boldsymbol{M}}|\times|\mathcal{I}_{\boldsymbol{M}}|}$,
then $\boldsymbol{T}$ is invertible and its inverse is given by

$\displaystyle\boldsymbol{T}^{-1}=\sum_{n=0}^{\infty}\left(\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}-\boldsymbol{T}\right)^{n}.$

Using this property for
$\boldsymbol{T}=\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A}$ we have

$\displaystyle\big{\|}(\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A})^{-1}\big{\|}_{2}=\left\|\sum_{n=0}^{\infty}\left(\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}-\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A}\right)^{n}\right\|_{2}=\left\|\sum_{n=0}^{\infty}\boldsymbol{\mathcal{E}}^{n}\right\|_{2}\leq\sum_{n=0}^{\infty}\left\|\boldsymbol{\mathcal{E}}^{n}\right\|_{2},$
(3.50)

in case that
$\big{\|}\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}-\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A}\big{\|}_{2}=\|\boldsymbol{\mathcal{E}}\|_{2}<1$.
Hence, by (3.49) and (3.50) we obtain

$\displaystyle\kappa_{2}(\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A})\leq\left(1+\|\boldsymbol{\mathcal{E}}\|_{2}\right)\cdot\left(\sum_{n=0}^{\infty}\left\|\boldsymbol{\mathcal{E}}^{n}\right\|_{2}\right).$
(3.51)

Additionally, we know that $|\varepsilon_{\boldsymbol{k}}|\leq\varepsilon$,
$\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{2M}}$, with some $\varepsilon>0$
and therefore

$\displaystyle\|\boldsymbol{\mathcal{E}}\|_{2}\leq\|\boldsymbol{\mathcal{E}}\|_{\mathrm{F}}=\sqrt{\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\,\sum_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}}}|\varepsilon_{\boldsymbol{\ell}-\boldsymbol{k}}|^{2}}\leq\sqrt{\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\,\sum_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}}}\varepsilon^{2}}=\varepsilon\,|\mathcal{I}_{\boldsymbol{M}}|.$
(3.52)

In other words, the correctness of (3.50) is ensured if
$\varepsilon\,|\mathcal{I}_{\boldsymbol{M}}|<1$. Since the spectral norm is a
sub-multiplicative norm, (3.52) also implies
$\|\boldsymbol{\mathcal{E}}^{n}\|_{2}\leq\|\boldsymbol{\mathcal{E}}\|_{2}^{n}\leq\left(\varepsilon\,|\mathcal{I}_{\boldsymbol{M}}|\right)^{n}$.
Consequently, we have

$\displaystyle\sum_{n=0}^{\infty}\left\|\boldsymbol{\mathcal{E}}^{n}\right\|_{2}\leq\sum_{n=0}^{\infty}\left(\varepsilon\,|\mathcal{I}_{\boldsymbol{M}}|\right)^{n}=\frac{1}{1-\varepsilon\,|\mathcal{I}_{\boldsymbol{M}}|}.$
(3.53)

Thus, combining (3.51), (3.52) and (3.53) yields the assertion (3.48). ∎

###  3.4 Connection to certain density compensation approaches from literature

In literature a variety of density compensation approaches can be found that
are concerned with the setting of bandlimited functions and make use of a
$\mathrm{sinc}$ transform

$\displaystyle\boldsymbol{C}\coloneqq\bigg{(}|\mathcal{I}_{\boldsymbol{M}}|\,\mathrm{sinc}\left(M\pi\left(\boldsymbol{x}_{j}-\boldsymbol{M}^{-1}\odot\boldsymbol{\ell}\right)\right)\bigg{)}_{j=1,\,\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}}}^{N}\
\in\mathbb{R}^{N\times|\mathcal{I}_{\boldsymbol{M}}|}$ (3.54)

instead of the Fourier transform (3.6). Namely, instead of directly using the
quadrature (3.23) for reconstruction, in these methods it is inserted into the
inverse Fourier transform (3.20), i. e.,

$\displaystyle f(\boldsymbol{x})$
$\displaystyle=\int\limits_{\left[-\frac{M}{2},\frac{M}{2}\right)^{d}}\hat{f}(\boldsymbol{v})\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{v}\boldsymbol{x}}\,\mathrm{d}\boldsymbol{v}\approx\sum_{j=1}^{N}w_{j}\,f(\boldsymbol{x}_{j})\int\limits_{\left[-\frac{M}{2},\frac{M}{2}\right)^{d}}\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{v}(\boldsymbol{x}_{j}-\boldsymbol{x})}\,\mathrm{d}\boldsymbol{v}$
$\displaystyle=\sum_{j=1}^{N}w_{j}\,f(\boldsymbol{x}_{j})\cdot|\mathcal{I}_{\boldsymbol{M}}|\,\mathrm{sinc}(M\pi(\boldsymbol{x}_{j}-\boldsymbol{x})),\quad\boldsymbol{x}\in\mathbb{R}^{d}.$
(3.55)

By using the $\mathrm{sinc}$ matrix
$\boldsymbol{C}\in\mathbb{R}^{N\times|\mathcal{I}_{\boldsymbol{M}}|}$ from
(3.54), the weight matrix $\boldsymbol{W}=\mathrm{diag}(w_{j})_{j=1}^{N}$ as
well as the vectors
$\boldsymbol{f}=\left(f(\boldsymbol{x}_{j})\right)_{j=1}^{N}$ and
$\boldsymbol{\tilde{f}}=(f(\boldsymbol{M}^{-1}\odot\boldsymbol{\ell}))_{\boldsymbol{l}\in\mathcal{I}_{\boldsymbol{M}}}$,
the evaluation of (3.4) at equispaced points
$\boldsymbol{M}^{-1}\odot\boldsymbol{\ell}$,
$\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}}$, can be denoted as
$\boldsymbol{\tilde{f}}\approx\boldsymbol{C}^{*}\boldsymbol{W}\boldsymbol{f}$.
Using the equispaced quadrature rule in (2.3), we find that evaluations
$\hat{f}(\boldsymbol{k})$ of (3.19) at the uniform grid points
$\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}$ can be approximated by (3.22)
by means of a simple FFT. In matrix-vector notation this can be written as
$\boldsymbol{\hat{f}}\approx\boldsymbol{\tilde{D}}^{*}\boldsymbol{F}_{|\mathcal{I}_{\boldsymbol{M}}|}^{*}\boldsymbol{\tilde{f}}$
where
$\boldsymbol{\hat{f}}=(\hat{f}(\boldsymbol{k}))_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}$,
$\boldsymbol{F}_{|\mathcal{I}_{\boldsymbol{M}}|}\coloneqq(\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}(\boldsymbol{M}^{-1}\odot\,\boldsymbol{\ell})})_{\boldsymbol{\ell},\,\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}$,
cf. (2.14), and
$\boldsymbol{\tilde{D}}\coloneqq\frac{1}{|\mathcal{I}_{\boldsymbol{M}}|}\,\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}$.
Thus, all in all one obtains an approximation of the form
$\boldsymbol{\hat{f}}\approx\boldsymbol{\tilde{D}}^{*}\boldsymbol{F}_{|\mathcal{I}_{\boldsymbol{M}}|}^{*}\boldsymbol{C}^{*}\boldsymbol{W}\boldsymbol{f}$.

Here some of these approaches, cf. [19], shall be reconsidered in the context
of the Fourier transform (3.6). We especially focus on the connection of the
approaches among each other as well as to our new method introduced in Section
3.1.

####  3.4.1 Density compensation using the pseudoinverse

Since (3.4) is in general not exactly solvable, we study the corresponding
least squares problem, instead, i. e., we look for the approximant that
minimizes the residual norm
$\big{\|}\boldsymbol{f}-\boldsymbol{A}\boldsymbol{\hat{f}}\big{\|}_{2}$. It is
known (e. g. [8, p. 15]) that this problem always has the unique solution

$\boldsymbol{\hat{f}}\approx\boldsymbol{\tilde{h}}^{\mathrm{pinv}}\coloneqq\boldsymbol{A}^{\dagger}\boldsymbol{f}$
(3.56)

with the Moore-Penrose pseudoinverse $\boldsymbol{A}^{\dagger}$. Comparing
(3.56) to the density compensation approach (3.6), the weights $w_{j}$ should
be chosen such that the matrix product $\boldsymbol{A}^{*}\boldsymbol{W}$
approximates the pseudoinverse $\boldsymbol{A}^{\dagger}$ as best as possible,
i. e., we study the optimization problem

$\underset{\boldsymbol{W}=\mathrm{diag}(w_{j})_{j=1}^{N}}{\text{Minimize }}\
\big{\|}\boldsymbol{A}^{*}\boldsymbol{W}-\boldsymbol{A}^{\dagger}\big{\|}_{\mathrm{F}}^{2},$
(3.57)

where $\|\cdot\|_{\mathrm{F}}$ denotes the Frobenius norm of a matrix. It was
shown in [62] that the solution to this least squares problem can be computed
as

$\displaystyle w_{j}$
$\displaystyle=\frac{[\boldsymbol{A}\boldsymbol{A}^{\dagger}]_{j,j}}{[\boldsymbol{A}\boldsymbol{A}^{*}]_{j,j}}=\frac{1}{|\mathcal{I}_{\boldsymbol{M}}|}\cdot[\boldsymbol{A}\boldsymbol{A}^{\dagger}]_{j,j},\quad
j=1,\dots N.$ (3.58)

However, since a singular value decomposition is necessary for the
calculations in (3.58), we obtain a high complexity of
$\mathcal{O}(N^{2}\,|\mathcal{I}_{\boldsymbol{M}}|+|\mathcal{I}_{\boldsymbol{M}}|^{3})$.
Therefore, we study some more sophisticated least squares approaches in the
following.

####  3.4.2 Density compensation using weighted normal equations of first kind

It is known, that every least squares solution to (3.4) satisfies the weighted
normal equations of first kind
$\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A}\boldsymbol{\hat{f}}=\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{f}$,
see e. g. [8, Thm. 1.1.2]. As already mentioned in Corollary 3.4, we have an
exact reconstruction formula (3.6) for all trigonometric polynomials (2.8) of
degree $\boldsymbol{M}$, if
$\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A}=\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}$
is fulfilled. Thus, we aim to compute optional weights $w_{j}$, $j=1,\dots,N$,
by considering the optimization problem

$\underset{\boldsymbol{W}=\mathrm{diag}(w_{j})_{j=1}^{N}}{\text{Minimize }}\
\big{\|}\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A}-\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}\big{\|}_{\mathrm{F}}^{2}.$
(3.59)

Analogous to [62] this could also be derived from (3.57) by introducing a
right-hand scaling in the domain of measured data and minimizing the Frobenius
norm of the weighted error matrix
$\boldsymbol{E}_{\mathrm{r}}\coloneqq\boldsymbol{E}\cdot\boldsymbol{A}$, where
$\boldsymbol{E}\coloneqq\boldsymbol{A}^{*}\boldsymbol{W}-\boldsymbol{A}^{\dagger}$
is the error matrix in (3.57).

In [59] it was shown that a solution
$\boldsymbol{W}=\operatorname*{diag}(\boldsymbol{w})$ to (3.59) can be
obtained by solving $\boldsymbol{S}\boldsymbol{w}=\boldsymbol{b}$ with

$\displaystyle\boldsymbol{S}$
$\displaystyle\coloneqq\Big{(}\left|\left[\boldsymbol{A}\boldsymbol{A}^{*}\right]_{j,h}\right|^{2}\Big{)}_{j,h=1}^{N}\quad\text{
and
}\quad\boldsymbol{b}=|\mathcal{I}_{\boldsymbol{M}}|\cdot\boldsymbol{1}_{N}.$
(3.60)

However, since $\boldsymbol{S}\boldsymbol{w}=\boldsymbol{b}$ is not separable
for single $w_{j}$, $j=1,\dots,N$, computing these weights is of complexity
$\mathcal{O}(N^{3})$. This is why the authors in [62] restricted themselves to
a maximal image size of $64\times 64$ pixels, which corresponds to setting
$M=64$.

####  3.4.3 Density compensation using weighted normal equations of second
kind

Another approach for density compensation factors is based on the weighted
normal equations of second kind

$\displaystyle\boldsymbol{A}\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{y}=\boldsymbol{f},\quad\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{y}=\boldsymbol{\hat{f}}.$
(3.61)

We recognize that by (3.61) we are given an exact approximation
$\boldsymbol{\hat{f}}=\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{f}$ of the
Fourier coefficients in (3.6) in case $\boldsymbol{y}=\boldsymbol{f}$, and
thereby $\boldsymbol{A}\boldsymbol{A}^{*}\boldsymbol{W}=\boldsymbol{I}_{N}$.
To this end, we consider the optimization problem

$\underset{\boldsymbol{W}=\mathrm{diag}(w_{j})_{j=1}^{N}}{\text{Minimize }}\
\|\boldsymbol{A}\boldsymbol{A}^{*}\boldsymbol{W}-\boldsymbol{I}_{N}\|_{\mathrm{F}}^{2}.$
(3.62)

As in Section 3.4.2, we remark that this optimization problem (3.62) could
also be derived from (3.57) by introducing an additional left-hand scaling in
the Fourier domain and minimizing the Frobenius norm of the weighted error
matrix
$\boldsymbol{E}_{\mathrm{l}}\coloneqq\boldsymbol{A}\cdot\boldsymbol{E}$.

######  Remark 3.16.

An analogous approach considering the $\mathrm{sinc}$ transform (3.54) instead
of the Fourier transform (3.6) was already studied in [52]. Another version
using a $\mathrm{sinc}$ transform evaluated at pointwise differences of the
nonequispaced points instead of (3.54) was studied in [12, 30], where it was
claimed that this approach coincides with the one in [52]. However, we remark
that due to the sampling theorem of Shannon-Whittaker-Kotelnikov, see (3.36),
applied to the function
$f(\boldsymbol{x})=\mathrm{sinc}(M\pi(\boldsymbol{x}_{j}-\boldsymbol{x}))$, i.
e.,

$\displaystyle\mathrm{sinc}(M\pi(\boldsymbol{x}_{j}-\boldsymbol{x}))=\sum_{\boldsymbol{\ell}\in\mathbb{Z}^{d}}\mathrm{sinc}(M\pi(\boldsymbol{x}_{j}-\boldsymbol{M}^{-1}\odot{\boldsymbol{\ell}}))\,\mathrm{sinc}\left(M\pi\left(\boldsymbol{x}-\boldsymbol{M}^{-1}\odot{\boldsymbol{\ell}}\right)\right)$

and its evaluation at $\boldsymbol{x}=\boldsymbol{x}_{h}$, $h=1,\dots,N$, this
claim only holds asymptotically for $|\mathcal{I}_{\boldsymbol{M}}|\to\infty$
in the setting of the $\mathrm{sinc}$ transform.

In contrast, when using the Fourier transform (3.6) this equality can directly
be seen. Then the analog to [30] utilizes an approximation of the form
$\boldsymbol{f}\approx\boldsymbol{H}\boldsymbol{W}\boldsymbol{f}$, where the
matrix $\boldsymbol{H}$ is defined as the system matrix of (3.4) evaluated at
pointwise differences of the nonequispaced points, i. e.,

$\boldsymbol{H}\coloneqq\Bigg{(}\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}(\boldsymbol{x}_{j}-\boldsymbol{x}_{h})}\Bigg{)}_{j,h=1}^{N}.$
(3.63)

Since by (2.7) we have $\boldsymbol{H}=\boldsymbol{A}\boldsymbol{A}^{*}$,
minimizing the approximation error

$\displaystyle\left\|\boldsymbol{H}^{*}\boldsymbol{W}\boldsymbol{f}-\boldsymbol{f}\right\|_{2}^{2}=\left\|\boldsymbol{A}\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{f}-\boldsymbol{f}\right\|_{2}^{2}$
$\displaystyle\leq\left\|\boldsymbol{A}\boldsymbol{A}^{*}\boldsymbol{W}-\boldsymbol{I}_{N}\right\|_{\mathrm{F}}^{2}\cdot\left\|\boldsymbol{f}\right\|_{2}^{2},$

leads to the optimization problem (3.62) as well.

It was shown in [52] that the minimizer of (3.62) is given by

$w_{j}=\frac{|\mathcal{I}_{\boldsymbol{M}}|}{\sum_{h=1}^{N}\big{|}\\!\left[\boldsymbol{A}\boldsymbol{A}^{*}\right]_{j,h}\\!\big{|}^{2}},\quad
j=1,\dots,N.$ (3.64)

Since for fixed $j$ the computation of
$\left[\boldsymbol{A}\boldsymbol{A}^{*}\right]_{j,h}$, $h=1,\dots,N,$ is of
complexity $\mathcal{O}(N\,|\mathcal{I}_{\boldsymbol{M}}|)$, the weights
(3.64) can be computed in $\mathcal{O}(N^{2}\,|\mathcal{I}_{\boldsymbol{M}}|)$
arithmetic operations. However, due to the explicit representation (3.63) the
computation of $\left[\boldsymbol{A}\boldsymbol{A}^{*}\right]_{j,h}$,
$h=1,\dots,N,$ for fixed $j$ can be accelerated by means of the NFFT (see
Algorithm 2.1). Then this step takes
$\mathcal{O}(|\mathcal{I}_{\boldsymbol{M}}|\log(|\mathcal{I}_{\boldsymbol{M}}|)+N)$
arithmetic operations and the overall complexity is given by
$\mathcal{O}(N\cdot|\mathcal{I}_{\boldsymbol{M}}|\log(|\mathcal{I}_{\boldsymbol{M}}|)+N^{2})$.

As mentioned in [52] one could also consider a simplified version of the
optimization problem (3.62) by reducing the number of conditions, e. g. by
summing the columns on both sides of
$\boldsymbol{A}\boldsymbol{A}^{*}\boldsymbol{W}=\boldsymbol{I}_{N}$ as

$\sum_{j=1}^{N}w_{j}\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}(\boldsymbol{x}_{h}-\boldsymbol{x}_{j})}=\sum_{j=1}^{N}\delta_{j,h}=1,\quad
h=1,\dots,N.$ (3.65)

By means of (2.7) this can be written as
$\boldsymbol{A}\boldsymbol{A}^{*}\boldsymbol{w}=\boldsymbol{1}_{N}$. Since
fast multiplication with $\boldsymbol{A}$ and $\boldsymbol{A}^{*}$ can be
realized using the NFFT (see Algorithm 2.1) and the adjoint NFFT (see
Algorithm 2.4), respectively, a solution to the linear system of equations
$\boldsymbol{A}\boldsymbol{A}^{*}\boldsymbol{w}=\boldsymbol{1}_{N}$ can be
computed iteratively with arithmetic complexity
$\mathcal{O}(|\mathcal{I}_{\boldsymbol{M}}|\log(|\mathcal{I}_{\boldsymbol{M}}|)+N)$.

Finally, we investigate the connection of this approach to our method
introduced in Section 3.1. To this end, suppose the linear system (3.12) is
fulfilled for given $\boldsymbol{w}\in\mathbb{C}^{N}$, i. e., by
$\boldsymbol{A}^{*}=\overline{\boldsymbol{A}^{T}}$ we have
$\left(\delta_{\boldsymbol{0},\boldsymbol{k}}\right)_{\boldsymbol{k}\in{\mathcal{I}_{\boldsymbol{M}}}}=\overline{\boldsymbol{A}^{T}\boldsymbol{w}}=\boldsymbol{A}^{*}\overline{\boldsymbol{w}}$.
Then multiplication with
$\boldsymbol{A}\in\mathbb{C}^{N\times|\mathcal{I}_{\boldsymbol{M}}|}$ in (2.7)
yields

$\boldsymbol{A}\boldsymbol{A}^{*}\overline{\boldsymbol{w}}=\boldsymbol{A}\cdot\left(\delta_{\boldsymbol{0},\boldsymbol{k}}\right)_{\boldsymbol{k}\in{\mathcal{I}_{\boldsymbol{M}}}}=\left(\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\delta_{\boldsymbol{0},\boldsymbol{k}}\cdot\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}}\right)_{j=1}^{N}=\boldsymbol{1}_{N}.$

In other words, an exact solution $\boldsymbol{w}$ to the linear system (3.12)
implies that the conjugate complex weights $\boldsymbol{\overline{w}}$ exactly
solve the system (3.65). However, the reversal does not hold true and
therefore (3.65) is not equivalent to (3.12). Moreover, we have seen in
Corollary 3.4 that an augmented variant of (3.12), namely (3.16), is necessary
to obtain an exact reconstruction
$\hat{f}_{\boldsymbol{k}}=h_{\boldsymbol{k}}^{\mathrm{w}}$ in (3.6) for
trigonometric polynomials (2.8) with maximum degree $\boldsymbol{M}$.

##  4 Direct inversion using matrix optimization

As seen in Remark 3.1, the previously considered density compensation
techniques can be regarded as an optimization of the sparse matrix
$\boldsymbol{B}\in\mathbb{R}^{N\times|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|}$
from the NFFT, cf. Section 2.1. Since density compensation allows only $N$
degrees of freedom, this limitation shall now be softened, i. e., instead of
searching for optimal scaling factors for the rows of $\boldsymbol{B}$, we now
study the optimization of each nonzero entry of the sparse matrix
$\boldsymbol{B}$, cf. [38]. To this end, we firstly have another look at the
equispaced setting. It is known by (3.1) and (3.2), that for equispaced points
and appropriately chosen parameters a one-sided inversion is given by
composition of the Fourier matrix and its adjoint. Hence, we aim to use this
result to find a good approximation of the inverse in the general setting.

Considering problem (3.4) we seek to find an appropriate matrix
$\boldsymbol{X}$ such that we have
$\boldsymbol{X}\boldsymbol{A}\approx\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}$,
since then we can simply compute an approximation of the Fourier coefficients
by means of
$\boldsymbol{X}\boldsymbol{f}=\boldsymbol{X}\boldsymbol{A}\boldsymbol{\hat{f}}\approx\boldsymbol{\hat{f}}$.
To find this left-inverse $\boldsymbol{X}$, we utilize the fact that in the
equispaced case it is known that (3.1) holds in the overdetermined setting
$|\mathcal{I}_{\boldsymbol{M}}|\leq N$. In addition, we also incorporate the
approximate factorization
$\boldsymbol{A}^{*}\approx\boldsymbol{D}^{*}\boldsymbol{F}^{*}\boldsymbol{B}^{*}$
of the adjoint NFFT, cf. Section 2.2, with the matrices
$\boldsymbol{D}\in\mathbb{C}^{|\mathcal{I}_{\boldsymbol{M}}|\times|\mathcal{I}_{\boldsymbol{M}}|}$,
$\boldsymbol{F}\in\mathbb{C}^{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|\times|\mathcal{I}_{\boldsymbol{M}}|}$
and
$\boldsymbol{B}\in\mathbb{R}^{N\times|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|}$
defined in (2.13), (2.14) and (2.15). Combining both ingredients we aim for an
approximation of the form
$\boldsymbol{D}^{*}\boldsymbol{F}^{*}\boldsymbol{B}^{*}\boldsymbol{A}\approx\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}$.
To achieve an approximation like this, we aim to modify the matrix
$\boldsymbol{B}$ such that its sparse structure with at most $(2m+1)^{d}$
entries per row and consequently the arithmetic complexity of its evaluation
is preserved. A matrix satisfying this property we call ${(2m+1)^{d}}$-sparse.

######  Remark 4.1.

We remark that this approach can also be deduced from the density compensation
method in Section 3 as follows. By Corollary 3.4 it is known that an exact
reconstruction needs to satisfy
$\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{A}=\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}$.
Since the reconstruction shall be realized efficiently by means of an adjoint
NFFT, one rather studies
$\boldsymbol{D}^{*}\boldsymbol{F}^{*}\boldsymbol{B}^{*}\boldsymbol{W}\boldsymbol{A}\approx\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}$.
Using the definition
$\boldsymbol{\tilde{B}}\coloneqq\boldsymbol{W}^{*}\boldsymbol{B}$ as in Remark
3.1, we end up with an approximation of the form
$\boldsymbol{D}^{*}\boldsymbol{F}^{*}\boldsymbol{\tilde{B}}^{*}\boldsymbol{A}\approx\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}$.
Thus, optimizing each nonzero entry of the sparse matrix $\boldsymbol{B}$
using this approximation is the natural generalization of the density
compensation method from Section 3.

Let $\boldsymbol{\tilde{B}}$ denote such a modified matrix. By defining
$\boldsymbol{\tilde{h}}:=\boldsymbol{D}^{*}\boldsymbol{F}^{*}\boldsymbol{\tilde{B}}^{*}\boldsymbol{f}$,
we recognize that the minimization of the approximation error

$\displaystyle\big{\|}\boldsymbol{\tilde{h}}-\boldsymbol{\hat{f}}\big{\|}_{2}$
$\displaystyle=\big{\|}\boldsymbol{D}^{*}\boldsymbol{F}^{*}\boldsymbol{\tilde{B}}^{*}\boldsymbol{f}-\boldsymbol{\hat{f}}\big{\|}_{2}=\big{\|}\boldsymbol{D}^{*}\boldsymbol{F}^{*}\boldsymbol{\tilde{B}}^{*}\boldsymbol{A}\boldsymbol{\hat{f}}-\boldsymbol{\hat{f}}\big{\|}_{2}$
$\displaystyle=\big{\|}\big{(}\boldsymbol{D}^{*}\boldsymbol{F}^{*}\boldsymbol{\tilde{B}}^{*}\boldsymbol{A}-\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}\big{)}\boldsymbol{\hat{f}}\big{\|}_{2}\leq\big{\|}\boldsymbol{D}^{*}\boldsymbol{F}^{*}\boldsymbol{\tilde{B}}^{*}\boldsymbol{A}-\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}\big{\|}_{\mathrm{F}}\,\big{\|}\boldsymbol{\hat{f}}\big{\|}_{2}$
(4.1)

implies the optimization problem

$\underset{\boldsymbol{\tilde{B}}\in\mathbb{R}^{N\times{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|}}\colon\boldsymbol{\tilde{B}}\,(2m+1)^{d}\text{-sparse
}}{\text{Minimize }}\
\big{\|}\boldsymbol{D}^{*}\boldsymbol{F}^{*}\boldsymbol{\tilde{B}}^{*}\boldsymbol{A}-\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}\big{\|}_{\mathrm{F}}^{2}.$
(4.2)

Note that a similar idea for the forward problem, i. e., the evaluation of
(2.5), was already studied in [51]. By the definition of the Frobenius norm we
have $\|\boldsymbol{Z}\|_{\mathrm{F}}=\|\boldsymbol{Z}^{*}\|_{\mathrm{F}}$,
such that (4.2) is equivalent to its adjoint

$\underset{\boldsymbol{\tilde{B}}\in\mathbb{R}^{N\times{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|}}\colon\boldsymbol{\tilde{B}}\,(2m+1)^{d}\text{-sparse
}}{\text{Minimize }}\
\big{\|}\boldsymbol{A}^{*}\boldsymbol{\tilde{B}}\boldsymbol{F}\boldsymbol{D}-\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}\big{\|}_{\mathrm{F}}^{2}.$
(4.3)

Since it is known by (2.14) that
$\boldsymbol{F}^{*}\boldsymbol{F}=|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|\,\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}$
and
$\boldsymbol{D}\in\mathbb{R}^{|\mathcal{I}_{\boldsymbol{M}}|\times|\mathcal{I}_{\boldsymbol{M}}|}$
is diagonal by (2.13), we have
$\tfrac{1}{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|}\boldsymbol{D}^{-1}\boldsymbol{F}^{*}\boldsymbol{F}\boldsymbol{D}=\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}$.
Thus, due to the fact that the Frobenius norm is a submultiplicative norm, we
have

$\displaystyle\big{\|}\boldsymbol{A}^{*}\boldsymbol{\tilde{B}}\boldsymbol{F}\boldsymbol{D}-\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}\big{\|}_{\mathrm{F}}$